function FSB(idat,scannum,reconum,fidopt,recox,recoy)

% FSB: Rapid preview fMRI toolbox
% GUI for performing various image manipulations to obtain better image
% quality from noisy fmri data.
%
% EXAMPLE:
% FSB('E06.FQ1',15,1,0)
% FSB(idat)
% FSB
%
% INPUT:
% idat          - image data/scan name
% scannum       - scan number
% reconum       - reconstruction number
% fidopt        - select the kind of image to be displayed
%               - If this parameter is set to 1 or 2, the raw image is read in, and
%               the reconstruction done in fMRI Sandbox
%               1: Display phase images
%               2: Display magnitude images
% recox         - reconstruction size in x if fidopt is one
% recoy         - reconstruction size in y
%
% OUTPUT:
% various output can be generated from within the GUI
%
% NOTES:
% Unfortunately, noisy fmri data have some properties that make them
% sometimes unsuitable for further analysis with standardized tools like SPM. With
% this toolbox, you may find the tools you need to prepare your data for
% further analysis. It allows you to see what kind of activation to expect from
% "real" statistics programs like SPM. For more precise analysis, please
% refer to other tools.
%
% Steffen Stoewer 05.09.08
% Graphics display and arrangement changed, new ROI display algorithm
% implemented.
% Please send comments and questions to stoewer@tuebingen.mpg.de
%
% Steffen Stoewer 05.01.10
% Completely restructured and modularized
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

warning off MATLAB:FINITE:obsoleteFunction
warning off MATLAB:stats:regress:RankDefDesignMat
warning off;
try
    feature('memstats')
catch
end
datascaling =1;

%~~~~~~~~~~~~~~~~~~~~~~
% set path to fMRI Sandbox directories - currently not necesssary
%~~~~~~~~~~~~~~~~~~~~~~
%addpath(genpath(cd));

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Enter here the path where your scanner files are located
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

global STDPATH acqp reco;
STDPATH.pv='W:\';

%~~~~~~~~~~~~~~~~~~~~~~~
% initialize variables
%~~~~~~~~~~~~~~~~~~~~~~~

mpimode = 0;
testmode =1;
maskit =1;
newmap = 0;

Map = struct;
sandbox = struct;
scaninfo = struct;
param = struct;

%~~~~~~~~~~~~~~~~~~~~~~~
% Create a struct to load files here
%~~~~~~~~~~~~~~~~~~~~~~~

if exist('idat','var')==1;
    loadstruct.idat = idat;
end

if exist('reconum','var')== 1;
    loadstruct.reconum = reconum;
end

if exist('scannum','var') ==1;
    loadstruct.scannum = scannum;
end

if exist('fidopt','var') == 1;
    loadstruct.fidopt = fidopt;
end

if exist('recox','var') == 1;
    loadstruct.recox = recox;
end

if exist('recoy','var') == 1;
    loadstruct.recoy = recoy;
end

if ~exist('loadstruct','var')
    loadstruct = struct;
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Load the config file with any parameters you deem important...
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

try
    [Map fsbdir mpimode testmode newmap maskit] =  sandbox_config;
    % sandbox_config;
    disp('Sandbox_config found...');
catch
    errmsg = lasterror;
    fsb_error_handler(errmsg);
    disp('No sandbox_config found, using generic parameters...');
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  Create and hide the GUI
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
sandbox_GUI = figure('Visible','off','Position',[0,0,990,720],'Resize','off','Toolbar','none','MenuBar','none','NumberTitle','off','Name','fMRI Sandbox');
movegui(sandbox_GUI,'northwest') % Upper left corner

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  GUI buttons, sliders and axes
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~
% File handling section:
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[20,500,120,215]);

uicontrol('Style','text',...
    'String','File','fontsize',10,...
    'ForegroundColor',[1 1 1],...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'Position',[30,690,100,17]);

uicontrol('Style','pushbutton','String','Load',...
    'Position',[30,660,100,25],...
    'TooltipString','Load multiple files at once',...
    'Callback',{@File_Load_Callback});

uicontrol('Style','pushbutton','String','Save',...
    'Position',[30,630,50,25],...
    'TooltipString','Save current data to disk',...
    'Callback',{@File_Save_Callback});

uicontrol('Style','pushbutton','String','Export',...
    'Position',[80,630,50,25],...
    'TooltipString','Save current data to disk',...
    'Callback',{@File_Export_Callback});

uicontrol('Style','pushbutton','String','Batch process',...
    'Position',[30,600,100,25],...
    'TooltipString','Run a batch process',...
    'Callback',{@Batch_processing_Callback});

uicontrol('Style','text','String','Load Design',...
    'Position',[30,580,100,12],...
    'TooltipString','Load design information');

uicontrol('Style','popupmenu',...
    'String',{'Auto','Blocked','Event related'},...
    'Position',[30,565,100,12],...
    'TooltipString','Select type of Design used',...
    'Callback',{@Design_load_Callback});


uicontrol('Style','text','String','SPM integration',...
    'Position',[30,540,100,12],...
    'TooltipString','Load design information');

uicontrol('Style','pushbutton',...
    'String','SPM selection GUI',...
    'Position',[30,515,100,25],...
    'TooltipString','Select SPM functions',...
    'Callback',{@SPM5_Button_Callback});

%~~~~~~~~~~~~~
% Artefact handling section
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[140,500,120,215]);

uicontrol('Style','text',...
    'String','Trial','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'ForegroundColor',[1 1 1],...
    'BackgroundColor',[0 0 1],...
    'Position',[150,690,100,17]);

uicontrol('Style','pushbutton','String','All',...
    'Position',[150,660,50,25],...
    'TooltipString','Select all trials',...
    'Callback',{@Trial_Sel_All_Callback});

uicontrol('Style','pushbutton','String','None',...
    'Position',[200,660,50,25],...
    'TooltipString','Deselect all Trials',...
    'Callback',{@Trial_Sel_None_Button_Callback});

uicontrol('Style','pushbutton','String','Remove Artefacts',...
    'Position',[150,642,100,13],...
    'TooltipString','Perform artefact removal',...
    'Callback',{@Artefact_Remove_Button_Callback});

h_detart_slider = uicontrol('Style','slider',...
    'Position',[150,630,100,12],...
    'Min',0,'Max',20,'Value',5,...
    'SliderStep',[(0.001) (0.1)],...
    'TooltipString','Preview and select artefact detection parameter',...
    'Callback',{@Artefact_detection_Slider_Callback});

uicontrol('Style','pushbutton','String','RemTrial',...
    'Position',[150,600,50,25],...
    'TooltipString','Remove selected trial periods',...
    'Callback',{@Trial_label_selected_Callback});

uicontrol('Style','pushbutton','String','LabVolume',...
    'Position',[200,600,50,25],...
    'TooltipString','Label selected volume as artefact',...
    'Callback',{@Volume_label_selected_Callback});

uicontrol('Style','pushbutton','String','Keep Selected',...
    'Position',[150,570,100,25],...
    'TooltipString','Keep selected trial(s)',...
    'Callback',{@Trial_Keep_Sel_Callback});

uicontrol('Style','pushbutton','String','Delete Selected',...
    'Position',[150,540,100,25],...
    'TooltipString','Delete selected trial(s)',...
    'Callback',{@Trial_Del_Sel_Callback});

uicontrol('Style','Pushbutton',...
    'String','Delete Volume',...
    'Position',[150,510,100,25],...
    'TooltipString','Delete selected volume',...
    'Callback',{@Volume_Delete_Selected_Callback});

%~~~~~~~~~~~~~
% Brain handling section
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[260,500,240,215]);

uicontrol('Style','text',...
    'String','Brain','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[270,690,220,17]);

% if mpimode ==1;
% 
%     uicontrol('Style','popupmenu',...
%         'String',{'Movement data','Toolbox','Unghost',...
%         'Realign 1D CoM','Realign 1D Corr','Realign 1D DFT',...
%         'Realign 2D Corr','Realign 2D Corr new','Realign Corr_AB',...
%         'Realign chisquare','Realign Resize',...
%         'Realign Demon 2D','Realign Demon 3D','Realign bspline 3D',...
%         'Realign Rigid Session','Realign Rigid full','Realign Rigid fine'},...
%         'Position',[270,660,100,25],...
%         'TooltipString','Realign volumes with various methods',...
%         'Callback',{@Volume_Realign_Callback});
% else

    uicontrol('Style','pushbutton',...
        'String','MotCor',...
        'Position',[270,660,50,25],...
        'TooltipString','Look for motion',...
        'Callback',{@Volume_Realign_Callback});
    
    uicontrol('Style','pushbutton',...
        'String','Toolbox',...
        'Position',[320,660,50,25],...
        'TooltipString','Realign datasets',...
        'Callback',{@Volume_Realign_Callback});

% end

uicontrol('Style','text','String','Apply Lens Filter',...
    'Position',[270,612,100,13],...
    'TooltipString','Apply a lens filter to facilitate brain extraction',...
    'Callback',{@Brain_lens_filter_Button_Callback});

h_lensfilt_slid = uicontrol('Style','Slider',...
    'Position',[270,600,100,12],...
    'Min',1,'Max',7,'Value',3,...
    'SliderStep',[(0.5/7) (1/7)],...
    'TooltipString','Select and preview lens filtering value',...
    'Callback',{@Brain_Lens_Filter_Callback});

uicontrol('Style','pushbutton','String','IntCor',...
    'Position',[270,630,50,25],...
    'TooltipString','Correct intensity over volume',...
    'Callback',{@Brain_Correct_Intensity_Gradient_Callback});

uicontrol('Style','pushbutton','String','Toolbox',...
    'Position',[320,630,50,25],...
    'TooltipString','Call Brain Extraction Toolbox',...
    'Callback',{@Brain_Extract_Toolbox_Callback});

uicontrol('Style','pushbutton','String','Extract',...
    'Position',[270,582,100,13],...
    'TooltipString','Perform brain extraction',...
    'Callback',{@Brain_Extract_Button_Callback});

h_extrbrain_slid = uicontrol('Style','slider',...
    'Position',[270,570,100,12],...
    'Min',1,'Max',12000,'Value',1,...
    'SliderStep',[(600/12000) (1200/12000)],...
    'TooltipString','Select and preview brain extraction parameter',...
    'Callback',{@Brain_Extract_Slider_Callback});

uicontrol('Style','pushbutton','String','Crop',...
    'Position',[270,540,100,25],...
    'TooltipString','Crop brain to reduce memory used and speed up processing',...
    'Callback',{@Brain_Crop_Callback});

uicontrol('Style','pushbutton','String','Pad',...
    'Position',[270,510,100,25],...
    'TooltipString','Pad brain to fit smallest common image size',...
    'Callback',{@Brain_Pad_Callback});

%~~~~~~~~~~~~~
% Normalization section
%~~~~~~~~~~~~~
if mpimode ==1

    uicontrol('Style','popupmenu',...
        'String',{'Process... ','Normalization Toolbox',...
        'Remove outliers','Normalize Volumes',...
        'Detrend scans', 'Detrend scans nonlinear','Interpolate Artefacts',...
        'Interpolate Artefacts and band-pass filter data',...
        'Interpolate Artefacts and high-pass filter data',...
        'Interpolate Artefacts and low-pass filter data',...
        'Enhance contrast',...
        'Normalize Sessions','Normalize Runs','Normalize Trials',...
        'Normalize Baseline multiplicative','Normalize Baseline divisive',...
        'Normalize Baseline additive','Normalize Baseline Z-score',...
        'Normalize Noise','Normalize ROI','Normalize Z-score'},...
        'Position',[390,660,100,25],...
        'TooltipString','Different normalization options',...
        'Callback',{@Volume_Normalize_Callback});
else

    uicontrol('Style','popupmenu',...
        'String',{'Process...',...
        'Remove outliers','Normalize Volumes',...
        'Detrend scans',...
        'Interpolate Artefacts and band-pass filter data',...
        'Interpolate Artefacts and high-pass filter data',...
        'Interpolate Artefacts and low-pass filter data',...
        'Normalize Sessions','Normalize Runs','Normalize Trials',...
        'Normalize Baseline additive',...
        'Normalize Noise','Normalize ROI','Normalize Z-score'},...
        'Position',[390,660,100,25],...
        'TooltipString','Temporally normalize data',...
        'Callback',{@Volume_Normalize_Callback});
end

uicontrol('Style','pushbutton',...
    'String','SNR Threshold',...
    'Position',[390,642,100,13],...
    'TooltipString','Threshold maps with SNR',...
    'Callback',{@Keep_Modifications});

uicontrol('Style','Slider',...
    'Position',[390,630,100,12],...
    'Min',10,'Max',50,'Value',10,...
    'SliderStep',[(0.1) (0.2)],...
    'TooltipString','Select SNR threshold',...
    'Callback',{@Threshold_SNR_Slider_Callback});

h_highres = uicontrol('Style','radiobutton',...
    'String','Anatomy Overlay',...
    'Position',[390,570,105,25],...
    'TooltipString','Use high resolution overlay image for map projections...',...
    'Callback',{@Overlay_use_Callback});

h_showslices = uicontrol('Style','radiobutton','String','Show slices',...
    'Position',[390,600,100,25],...
    'TooltipString','Slice display on/off',...
    'Callback',{@Slice_display_Radiobutton_Callback});

uicontrol('Style','pushbutton','String','Smooth Temporal',...
    'Position',[390,552,100,13],...
    'TooltipString','Temporally smooth data',...
    'Callback',{@Smooth_Idat_Temporal_Button_Callback});

h_tsmooth_slider = uicontrol('Style','Slider',...
    'Position',[390,540,100,12],...
    'Min',1,'Max',7,'Value',2,...
    'SliderStep',[(1/6) (2/6)],...
    'TooltipString','Select temporal smoothing kernel',...
    'Callback',{@Smooth_Idat_Temporal_Slider_Callback});

uicontrol('Style','pushbutton','String','Smooth Spatial',...
    'Position',[390,522,100,13],...
    'TooltipString','Perform spatial smoothing',...
    'Callback',{@Filter_Idat_Spatial_Callback});

h_ssmooth_slider = uicontrol('Style','Slider',...
    'Position',[390,510,100,12],...
    'Min',0,'Max',10,'Value',2,...
    'SliderStep',[(0.1) (0.2)],...
    'TooltipString','Select spatial smoothing kernel',...
    'Callback',{@Filter_Idat_Spatial_Slider_Callback});

%~~~~~~~~~~~~~
% Map section
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[500,500,240,215]);

uicontrol('Style','text',...
    'String','Map','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[510,690,210,17]);

uicontrol('Style','popupmenu',...
    'String',{'Correlation','GLM','Robust Regression','ICA'},...
    'Position',[510,600,100,25],...
    'TooltipString','Select statistical model',...
    'Callback',{@Map_model_selection_Callback});

h_calcorrelation2_slid = uicontrol('Style','Slider',...
    'Position',[510,570,100,12],...
    'Min',0,'Max',1,'Value',0.2,...
    'SliderStep',[(0.001/1) (0.1/1)],...
    'TooltipString','Select correlation threshold',...
    'Callback',{@Map_threshold_slider_Callback});

uicontrol('Style','text','String','Map Threshold',...
    'Position',[510,582,100,13]);

uicontrol('Style','text','String','Correlation : ',...
    'Position',[510,552,60,13]);

h_probability_button = uicontrol('Style','edit','String','0.0000',...
    'Position',[570,552,40,13],...
    'TooltipString','Select correlation threshold',...
    'Callback',{@Map_threshold_edit_Callback});

h_probability_slid = uicontrol('Style','text','String','Probability : 0.0000',...
    'Position',[510,540,100,12]);

uicontrol('Style','text','String','Cluster threshold',...
    'Position',[510,522,100,13]);

h_clustering_slid = uicontrol('Style','Slider',...
    'Position',[510,510,100,12],...
    'Min',1,'Max',101,'Value',1,...
    'SliderStep',[(10/100) (20/100)],...
    'TooltipString','Select minimum cluster size',...
    'Callback',{@Map_cluster_size_Slider_Callback});

uicontrol('Style','text','String','Predictor',...
    'Position',[510,635,50,15]);

predsel = uicontrol('Style','edit',...
    'String','1',...
    'Position',[560,630,50,25],...
    'TooltipString','Enter the predictor you want to display',...
    'Callback',{@Map_predictor_selection_Callback});

uicontrol('Style','pushbutton','String','Calculate full map',...
    'Position',[510,660,100,25],...
    'TooltipString','Calculate a statistical map for whole volume',...
    'Callback',{@Map_calculate_volume_Callback});

h_mapdisp = uicontrol('Style','radiobutton','String','Map display on',...
    'Position',[620,660,100,25],...
    'TooltipString','Display statistical map for predictor chosen',...
    'Callback',{@Map_display_Button_Callback});

if mpimode == 1

    uicontrol('Style','popupmenu',...
        'String',{'Convert trial structure','Interpolate artefacts and trial data','Create regressor for current ROI'},...
        'Position',[620,630,100,25],...
        'TooltipString','Process your dataset for baseline analysis',...
        'Callback',{@Baseline_conversion_Callback});
end

h_trialselrb = uicontrol('Style','radiobutton','String','Only selected',...
    'Position',[620,600,100,25],...
    'TooltipString','Analyze only selected trial(s)',...
    'Callback',{@Trial_analyze_selected_Callback});

h_trialselrb = uicontrol('Style','radiobutton','String','Only selected',...
    'Position',[620,600,100,25],...
    'TooltipString','Analyze only selected trial(s)',...
    'Callback',{@Trial_analyze_selected_Callback});

h_overlaytrial = uicontrol('Style','radiobutton','String','Trial Timecourse',...
    'Position',[620,570,100,25],...
    'TooltipString','Trial-locked hemodynamic timecourse',...
    'Callback',{@Trial_display_timecourse_Button_Callback});

h_overlayclus = uicontrol('Style','radiobutton','String','Cluster Data',...
    'Position',[620,540,100,25],...
    'TooltipString','Select Cluster level data display',...
    'Callback',{@Trial_TC_voxel_cluster_Callback});

uicontrol('Style','text','String','Roi Size',...
    'Position',[620,522,100,13]);

h_roi_slid = uicontrol('Style','slider',...
    'Position',[620,510,100,12],...
    'Min',1,'Max',100,'Value',1,...
    'SliderStep',[(1/100) (10/100)],...
    'TooltipString','Select and preview ROI size',...
    'Callback',{@ROI_size_Slider_Callback});

%~~~~~~~~~~~~~
% Exploratory Section
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[730,500,120,215]);

uicontrol('Style','text',...
    'String','Explore','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[740,690,100,17]);

uicontrol('Style','pushbutton','String','Reset HRF Shift',...
    'Position',[740,672,100,13],...
    'TooltipString','Reset the HRF shift to zero',...
    'Callback',{@HRF_shift_Button_Callback});

h_shifthrfslider = uicontrol('Style','slider',...
    'Position',[740,660,70,12],...
    'Min',-2,'Max',+2,'Value',0,...
    'SliderStep',[(0.1/4) (0.5/4)],...%'SliderStep',[(2/2) (2/2)],...
    'TooltipString','Select the amount of hrf shifting',...
    'Callback',{@HRF_shift_Slider_Callback});

h_shifthrfedit = uicontrol('Style','edit',...
    'String','0',...
    'Position',[810,660,30,12],...
    'TooltipString','Select the amount of hrf shifting',...
    'Callback',{@HRF_edit_Callback});

uicontrol('Style','pushbutton','String','Smooth HRF',...
    'Position',[740,642,100,13],...
    'TooltipString','Smooth the HRF model',...
    'Callback',{@HRF_smooth_Button_Callback});

h_smoothhrfslider = uicontrol('Style','slider',...
    'Position',[740,630,100,12],...
    'Min',0,'Max',5,'Value',0,...
    'SliderStep',[(1/4) (2/4)],...
    'TooltipString','Select the amount of smoothing',...
    'Callback',{@HRF_smooth_Slider_Callback});

uicontrol('Style','Pushbutton',...
    'String','AddBase',...
    'Position',[740,600,50,25],...
    'TooltipString','Add a Baseline volume to every trial',...
    'Callback',{@Trial_add_baseline_Callback});

uicontrol('Style','Pushbutton',...
    'String','RemBase',...
    'Position',[790,600,50,25],...
    'TooltipString','Remove a baseline volume from every trial',...
    'Callback',{@Trial_rem_baseline_Callback});

uicontrol('Style','Pushbutton',...
    'String','Fine-Tune',...
    'Position',[740,570,100,25],...
    'TooltipString','Fine-Tune regressor fit for ROI...',...
    'Callback',{@HRF_auto_adjust_Callback});

if mpimode ==1;
    uicontrol('Style','popupmenu','String',{'Spectrum','High-Pass Filter',...
        'Low-Pass Filter','Band-Pass Filter','Band-Pass Filter run by run',...
        'Remove discontinuities',...
        'Interpolate Artefacts and band-pass filter data',...
        'Interpolate Artefacts and low-pass filter data',...
        'Interpolate Artefacts and high-pass filter data',...
        'Interpolate out of trial periods and filter data'},...
        'Position',[740,540,100,25],...
        'TooltipString','Temporally filter image data',...
        'Callback',{@Filter_Temporal_Callback});
end

uicontrol('Style','text',...
    'String','High(s)', 'Position',[745,525,50,15]);

h_filter_high = uicontrol('Style','edit',...
    'Position',[740,505,50,20],...
    'String','40',...
    'TooltipString','Edit High Pass Filter Frequency',...
    'Callback',{@Filter_high_Frequency_Callback});

uicontrol('Style','text',...
    'String','Low(s)', 'Position',[785,525,50,15]);

h_filter_low = uicontrol('Style','edit',...
    'Position',[790,505,50,20],...
    'String','5',...
    'TooltipString','Edit Low Pass Filter Frequency',...
    'Callback',{@Filter_low_Frequency_Callback});

%~~~~~~~~~~~~~
% General tabs
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[860,420,120,350]);

uicontrol('Style','text',...
    'String','General','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[870,690,100,17]);

uicontrol('Style','pushbutton','String','Help',...
    'fontweight','b','horizontalalignment','center',...
    'Position',[870,660,100,25],'BackgroundColor',[1 0 0],...
    'ForegroundColor',[1 1 1],...
    'TooltipString','Read the helpfile',...
    'Callback',{@Help_Button_Callback});

uicontrol('Style','pushbutton','String','Undo',...
    'Position',[870,630,100,25],...
    'TooltipString','Undo last operation',...
    'Callback',{@Undo_Button_Callback});

uicontrol('Style','text','String','Color',...
    'Position',[870,605,50,15]);

uicontrol('Style','popupmenu',...
    'String',{'Gray','Hot','Jet','HSV','Cool','Lines','Copper','Flag','Prism','Bone'},...
    'Position',[920,600,50,25],...
    'TooltipString','Select a color scheme for display',...
    'Callback',{@Colormap_Selection_Callback});

%~~~~~~~~~~~~~
% Write Operations
%~~~~~~~~~~~~~

uicontrol('Style','Pushbutton',...
    'String','Write Mean Brain',...
    'Position',[870,575,100,25],...
    'TooltipString','Write mean brain...',...
    'Callback',{@Volume_Write_Mean_Callback});

uicontrol('Style','Pushbutton',...
    'String','Write Selected',...
    'Position',[870,545,100,25],...
    'TooltipString','Write selected volume...',...
    'Callback',{@Volume_Write_Selected_Callback});

uicontrol('Style','text',...
    'String','FlipVol',...
    'Position',[870,510,50,25]);

uicontrol('Style','popupmenu',...
    'String',{'no','x','y','z','rotate right','rotate left','shiftdim1','shiftdim2','shiftdim3','Keep transform'},...
    'Position',[920,515,50,25],...
    'TooltipString','Flip the volume displayed',...
    'Callback',{@Volume_flip_selection_Callback});

uicontrol('Style','Pushbutton',...
    'String','Del ROI',...
    'Position',[870,490,50,25],...
    'TooltipString','Overwrite ROI area with zeros...',...
    'Callback',{@ROI_voxels_delete_Callback});

uicontrol('Style','Pushbutton',...
    'String','Keep ROI',...
    'Position',[920,490,50,25],...
    'TooltipString','Keep only ROI area...',...
    'Callback',{@ROI_voxels_keep_Callback});

uicontrol('Style','Pushbutton',...
    'String','Write TC',...
    'Position',[870,460,100,25],...
    'TooltipString','Write ROI/Cluster time course...',...
    'Callback',{@TC_write_Cluster_Callback});

uicontrol('Style','popupmenu',...
    'String',{'int8','int16','single','double'},...
    'Position',[870,430,50,25],...
    'TooltipString','Change datatype of image data...',...
    'Callback',{@datatype_Callback});

h_transform_display = uicontrol('Style','radiobutton',...
    'String','MAT',...
    'Position',[920,430,50,25],...
    'TooltipString','Display transforms...',...
    'Callback',{@transform_display_Callback});

%~~~~~~~~~~~~~~~~~~~~~~~~~
% Information
%~~~~~~~~~~~~~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[860,310,120,105]);

uicontrol('Style','text',...
    'String','Information','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[870,390,100,17]);

uicontrol('Style','text',...
    'String','xyz= ',...
    'Position',[865,355,25,25]);

h_xhdr = uicontrol('Style','edit',...
    'Position',[890,360,25,25],...
    'String','1',...
    'TooltipString','Edit Voxelsize x',...
    'Callback',{@HDR_Callback});

h_yhdr = uicontrol('Style','edit',...
    'Position',[918,360,25,25],...
    'String','1',...
    'TooltipString','Edit Voxelsize y',...
    'Callback',{@HDR_Callback});

h_zhdr = uicontrol('Style','edit',...
    'Position',[946,360,25,25],...
    'String','1',...
    'TooltipString','Edit Voxelsize z',...
    'Callback',{@HDR_Callback});

uicontrol('Style','text',...
    'String','orig=',...
    'Position',[865,325,25,25]);

h_xorig = uicontrol('Style','edit',...
    'Position',[890,330,25,25],...
    'String','1',...
    'TooltipString','Edit Origin x',...
    'Callback',{@HDR_Callback});

h_yorig = uicontrol('Style','edit',...
    'Position',[918,330,25,25],...
    'String','1',...
    'TooltipString','Edit Origin y',...
    'Callback',{@HDR_Callback});

h_zorig = uicontrol('Style','edit',...
    'Position',[946,330,25,25],...
    'String','1',...
    'TooltipString','Edit Origin z',...
    'Callback',{@HDR_Callback});

%~~~~~~~~~~~~~
% Diagnostics
%~~~~~~~~~~~~~

uicontrol('Style','frame',...
    'Position',[860,180,120,125]);

uicontrol('Style','text',...
    'String','Diagnostics','fontsize',10,...
    'fontweight','b','horizontalalignment','center',...
    'BackgroundColor',[0 0 1],...
    'ForegroundColor',[1 1 1],...
    'Position',[870,280,100,17]);

uicontrol('Style','popupmenu',...
    'String',{'Scaninfo','Plot slices','Voxel Time Course','SNR Map',...
    'ROI SNR Map','CV Map','TN Map','STD Map','Z-score Map','Special Map',...
    'AB Map','Ghost Map','Noise Time Course','Sandmap','STD Map2','Volume Histogram','Euclidian Distance'},...
    'Position',[870,250,100,25],...
    'TooltipString','Select diagnostic method here',...
    'Callback',{@Slice_Display_Callback});

if testmode ==1
    uicontrol('Style','Pushbutton',...
        'String','ICA',......
        'Position',[870,220,50,25],...
        'TooltipString','Run 10 component ICA over data',...
        'Callback',{@ICA_Callback});
    uicontrol('Style','Pushbutton',...
        'String','Toolbox',......
        'Position',[920,220,50,25],...
        'TooltipString','Call fastICA toolbox',...
        'Callback',{@ICA_Toolbox_Callback});

    uicontrol('Style','Pushbutton',...
        'String','Compress slices',......
        'Position',[870,190,100,25],...
        'TooltipString','Compress your slices 2>1',...
        'Callback',{@compress_slices_Callback});
end
%~~~~~~~~~~~~~
% Slider controls
%~~~~~~~~~~~~~

h_xaxis = uicontrol('Style','slider','String','x-axis',...
    'Position',[30,374,262,12],...
    'Min',1,'Max',100,'Value',50,...
    'SliderStep',[(1/100) (10/100)],...
    'Callback',{@Slider_x_Callback});

h_yaxis = uicontrol('Style','slider','String','y-axis',...
    'Position',[30,254,262,12],...
    'Min',1,'Max',100,'Value',50,...
    'SliderStep',[(1/100) (10/100)],...
    'Callback',{@Slider_y_Callback});

h_zaxis = uicontrol('Style','slider','String','z-axis',...
    'Position',[580,269,12,222],...
    'Min',1,'Max',100,'Value',50,...
    'SliderStep',[(1/100) (10/100)],...
    'Callback',{@Slider_z_Callback});

h_taxis = uicontrol('Style','Slider',...
    'String','t-axis',...
    'Position',[30,10,850,12],...
    'Min',1,'Max',100,'Value',50,...
    'SliderStep',[(1/100) (10/100)],...
    'Callback',{@Slider_t_Callback});

%~~~~~~~~~~~~~
% Define subplots
%~~~~~~~~~~~~~

hc = axes('Units','Pixels','DataAspectRatio',[1 1 1],'Position',[320,270,260,220]);axis off;
hb = axes('Units','Pixels','DataAspectRatio',[1 1 1],'Position',[30,270,260,100]);axis off;
xlabel(' ');
ha = axes('Units','Pixels','DataAspectRatio',[1 1 1],'Position',[30,390,260,100]);axis off;
colmap = axes('Units','Pixels','DataAspectRatio',[1 1 1],'Position',[312,270,5,220],'xTickLabel',' ');axis off;
hd = axes('Units','Pixels','DataAspectRatio',[1 1 1],'Position',[30,40,820,170]);
he = axes('Units','Pixels','Position',[630,270,220,210],'Color',[0.9 0.9 0.9]);  axis off;
title('% Signal change');xlabel('Trial timepoint');
txt_win = axes('Units','Pixels','Position',[870,120,150,160]); axis off;
beta_win = axes('Units','Pixels','Position',[870,20,120,100]); axis off;
info_win = axes('Units','Pixels','Position',[870,300,50,120]); axis off;
feedb_win = axes('Units','Pixels','Position',[30,240,200,25]); axis off;

%~~~~~~~~~~~~~~~~~
% Load file if asked to
%~~~~~~~~~~~~~~~~~

try
    [idat sandbox scaninfo] = fsb_preload_sandbox(loadstruct);
catch
    disp('****************************');
    disp('File not found');
    disp('****************************');
    return
end

clear loadstruct;

set(h_xhdr,'String',num2str(sandbox.hdr.dime.pixdim(2)));
set(h_yhdr,'String',num2str(sandbox.hdr.dime.pixdim(3)));
set(h_zhdr,'String',num2str(sandbox.hdr.dime.pixdim(4)));

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Determine size of the image
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

dim = size(idat);
Tempvar.slice_n = round(dim/2);
Tempvar.slice_n(4) = 1;

if size(dim,2)>3 && dim(4)==1;
    set(hd,'visible','off');
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Declare additional variables
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%Parameters
Tempvar.threshold_SNR = 10;
Tempvar.extr_param = 1;
Tempvar.detart_param = 3;
Tempvar.t_param = 1;

Tempvar.ssmooth_vox = 2;
Tempvar.ssmooth_param = Tempvar.ssmooth_vox/2.355;
Tempvar.roi =1;
Map.corparam = 0.2;
Map.probparam = 0.1;
Tempvar.feedb = 'Welcome to fMRI Sandbox';
Tempvar.do = 1;
Map.statsel = 'corrc';
Tempvar.clustparam =1;
Tempvar.lensfilt_param =1;
Tempvar.imagechange = 1;
Tempvar.slice_l = Tempvar.slice_n;
Tempvar.slice_h = Tempvar.slice_n(1:3);
Tempvar.dimratio = [1 1 1];

% Variables with at value of zero
Tempvar.overlayt = 0;
Tempvar.hoverlay = 0;
Tempvar.loverlay =0;
Tempvar.sparecor =0;
Tempvar.snr = 0;
Tempvar.mapdisp =0;
Tempvar.disp_disc = 0;
Tempvar.shiftparam = 0;
Tempvar.smoothed=0 ;
Tempvar.smoothparam=0;
Tempvar.ctr_disp2 =0;
Tempvar.intrial_nr =0;
Map.rc =0;
param.dispslices = 0;
Tempvar.psc =0;
Tempvar.mapcalc = 0;

if exist('hdr','var') ==0
    hdr = {};
else
    sandbox.hdr.dime.pixdim(2) = hdr.dime.pixdim(2);
    sandbox.hdr.dime.pixdim(3)=  hdr.dime.pixdim(3);
    sandbox.hdr.dime.pixdim(4) = hdr.dime.pixdim(4);
end

if exist('sandbox','var')==0 ||~isfield(sandbox,'hdr')
    sandbox.hdr.dime.pixdim(2) = 1;
    sandbox.hdr.dime.pixdim(3)=  1;
    sandbox.hdr.dime.pixdim(4) = 1;
end

hdr.dime.pixdim = [0 sandbox.hdr.dime.pixdim(2) sandbox.hdr.dime.pixdim(3) sandbox.hdr.dime.pixdim(4) 0 0 0 0];

%~~~~~~~~~~~~~~~~
% empty variables
%~~~~~~~~~~~~~~~~
otto =[];trial_hrf2_normalized=[];
trial_normalized = [];
hem2trial_normalized = [];
best=[]; bint=[];
hrf_disp =[]; hrf_disp2 = 0; param.hrf_pred =1;
analtrial=0; clustselect=0;
idat_red = []; hemodynamics_red =[]; temp =[];
r_cor_red = zeros(dim(1),dim(2),dim(3),1); r_cor_neg =[];
bhemodynamics =[]; hemodynamics2 =[];
hemodynamics_deriv =[]; dummy = [];
disp_onset=1;
trial_hrf=[];
trial_idat=[];
ctr_disp=[]; ctr_mass =[];
dctr_mass_xyz =[]; dctr_mass_r =[];
hidat =[]; hidat_spare = [];
hrfcalc.x = zeros(1,dim(1));
hrfcalc.y = zeros(1,dim(2));
hrfcalc.z = zeros(1,dim(3));

%~~~~~~~~~~~~~~~~
% Calculated parameters
%~~~~~~~~~~~~~~~~

hrf_vol_disp = squeeze(mean(mean(mean(idat(1:end,1:end,1:end,:)))));
min_hrf_vol = min(hrf_vol_disp);
max_hrf_vol = max(hrf_vol_disp);
avg_hrf_vol = mean(hrf_vol_disp);
psc_hrf_vol = (max_hrf_vol-min_hrf_vol)/avg_hrf_vol*100;
psc_vol_string = sprintf('%1.2f',psc_hrf_vol);
imean = fsb_calculate_mean(idat); % calculate average brain

%~~~~~~~~~~~~~~~~
% calculate image average value to use as threshold
%~~~~~~~~~~~~~~~~
imean_avg = mean(mean(mean(imean(round(Tempvar.slice_n(1)-Tempvar.slice_n(1)/6):round(Tempvar.slice_n(1)+Tempvar.slice_n(1)/6),...
    round(Tempvar.slice_n(2)-Tempvar.slice_n(2)/6):round(Tempvar.slice_n(2)+Tempvar.slice_n(2)/6),...
    round(Tempvar.slice_n(3)-Tempvar.slice_n(3)/6):round(Tempvar.slice_n(3)+Tempvar.slice_n(3)/6)))));

%conditional variables

if ~isfield(sandbox,'hemodynamics');
    sandbox.hemodynamics = [];
    sandbox.stats.r_cor = zeros(dim(1),dim(2),dim(3),1);
    disp_onset = 0;

else

    sandbox.stats.r_cor = zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
    imean2 = repmat(imean,[1 1 1 size(sandbox.hemodynamics,2)]);
    sandbox.stats.r_cor(imean2<imean_avg/3)=NaN; % Mask sandbox.stats.r_cor and all other map matrices with a threshold for brain/nobrain
end

r_cor_red = sandbox.stats.r_cor;
r_cor_neg = sandbox.stats.r_cor;
sandbox.stats.brob = sandbox.stats.r_cor;
sandbox.stats.b_est = sandbox.stats.r_cor;
sandbox.stats.p_rob= sandbox.stats.r_cor;
sandbox.stats.rob_est= sandbox.stats.r_cor;
sandbox.stats.rob_prob= sandbox.stats.r_cor;
sandbox.stats.glm_prob= sandbox.stats.r_cor;
sandbox.stats.t_test= sandbox.stats.r_cor;
hemodynamics_spare = sandbox.hemodynamics;
r_cor_spare =sandbox.stats.r_cor;
r_cor_red_spare = r_cor_red;
r_cor_neg_spare =r_cor_neg;
sandbox.stats.r_cor(idat(:,:,:,1) == 0) = NaN;

if ~isfield(sandbox,'intrial')
    if exist('intrial','var')==0;
        sandbox.intrial = [];
    else
        sandbox.intrial = intrial;
    end
end

if ~isempty(sandbox.intrial)
    if size(sandbox.intrial,2)<3;
        intrial2 = sandbox.intrial;
        intrial2(intrial2==0.01)=0;
        sandbox.intrial(:,3) = bwlabeln(intrial2(:,1));
    end
end

if ~isfield(sandbox,'onset')
    sandbox.onset = [];
end

intrial_spare = sandbox.intrial;
onset_spare = sandbox.onset;

b_idat = idat(:,:,:,1);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  Get the information to be displayed
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

slicez_data = squeeze(idat(:,:,Tempvar.slice_n(3),Tempvar.slice_n(4)));
slicez_data = slicez_data';
slicex_data = squeeze(idat(Tempvar.slice_n(1),:,:,Tempvar.slice_n(4)));
slicex_data = slicex_data';
slicey_data = squeeze(idat(:,Tempvar.slice_n(2),:,Tempvar.slice_n(4)));
slicey_data = slicey_data';
slicet_data = squeeze(idat(:,Tempvar.slice_n(2),Tempvar.slice_n(3),:));
idat_spare = idat;
refresh_all;
colormap gray;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Initialize the GUI
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

set(sandbox_GUI,'Visible','on');
refresh_figure;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
% Colormap menu
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function Colormap_Selection_Callback(source,eventdata)

        str = get(source, 'String');
        val = get(source,'Value');
        colormap(str{val});
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Predictor menu
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Map_predictor_selection_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        param.hrf_pred = get(predsel, 'String');
        param.hrf_pred = str2num(param.hrf_pred);

        % if more than one predictor was entered
        if size(param.hrf_pred,2) >1
            sandbox.hemodynamics(:,26)=sandbox.hemodynamics(:,param.hrf_pred(1))+sandbox.hemodynamics(:,param.hrf_pred(2));
            sandbox.onset(:,26)=sandbox.onset(:,param.hrf_pred(1))+sandbox.onset(:,param.hrf_pred(2));
            param.hrf_pred = 26;
            set(predsel,'String',num2str(26));
            Tempvar.feedb = ('Combined predictor: 66');
        end

        if param.hrf_pred>size(sandbox.hemodynamics,2)
            if param.hrf_pred==99
                %param.hrf_pred =1;
                param.hrf_pred = 25;
                sandbox.hemodynamics(:,25) = Tempvar.roimean;

                set(predsel,'String',num2str(param.hrf_pred));
                disp ('defaulting to ROI predictor');
                Tempvar.feedb = ('defaulting to ROI predictor');
            else
                param.hrf_pred=1;
                set(predsel,'String',num2str(param.hrf_pred));
                disp ('defaulting to first predictor');
                Tempvar.feedb = ('defaulting to first predictor');
            end

        else
            set(predsel,'String',num2str(param.hrf_pred));
            Tempvar.feedb = (['selected predictor: ' num2str(param.hrf_pred)]);
        end

        refresh_all;
        refresh_feedb;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
%Statistical model selection
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Map_model_selection_Callback(source,eventdata)

        str = get(source, 'String');
        val = get(source,'Value');

        switch str{val};

            case 'GLM'
                Map.statsel='glm';
            case 'ICA'
                Map.statsel='ICA';
            case 'Correlation'
                Map.statsel='corrc';
            case 'Robust Regression';
                Map.statsel='robreg';
            case 'T-Test';
                Map.statsel ='ttest';
            case 'ConMan';
                Map.statsel ='conman';

        end

        Tempvar.mapcalc = 0;
        Map_threshold_slider_Callback;
        refresh_all;

    end

    function Volume_flip_selection_Callback(source,eventdata)
        str = get(source, 'String');
        val = get(source,'Value');
        hidat = hidat_spare;

        switch str{val};
            case 'no'
                idat = idat_spare;
            case 'x'
                idat = flipdim(idat_spare,1);
                hidat = flipdim(hidat_spare,1);
            case 'y'
                idat = flipdim(idat_spare,2);
                hidat = flipdim(hidat_spare,2);
            case 'z'
                idat = flipdim(idat_spare,3);
                hidat = flipdim(hidat_spare,3);
            case 'rotate right'
                idat = permute(idat_spare,[2 1 3 4]);
                idat = flipdim(idat,1);
                hidat = permute(hidat_spare,[2 1 3 4]);
                hidat = flipdim(hidat,1);
            case 'rotate left'
                idat = permute(idat_spare,[2 1 3 4]);
                hidat = permute(hidat_spare,[2 1 3 4]);
            case 'shiftdim1'
                idat = permute(idat_spare,[2 1 3 4]);
                hidat = permute(hidat_spare,[2 1 3 4]);
            case 'shiftdim2'
                idat = permute(idat_spare,[3 2 1 4]);
                hidat = permute(hidat_spare,[2 1 3 4]);
            case 'shiftdim3'
                idat = permute(idat_spare,[2 3 1 4]);
                hidat = permute(hidat_spare,[2 3 1 4]);
            case 'Keep transform'
                Keep_Modifications;
                refresh_all;

        end

        dim = size(idat);
        dim2 = size(hidat);
        Tempvar.slice_n = [round(dim(1:3)/2) 1];
        Tempvar.slice_l = Tempvar.slice_n;
        if Tempvar.hoverlay ==1;
            Tempvar.slice_h = [round(dim2(1:3)/2) 1];
        end
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Crosshair buttondown callbacks
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Click_Time_Callback(source,eventdata)

        xy = get(gca,'CurrentPoint');
        sel_typ = get(gcbf,'SelectionType');
        switch sel_typ
            case 'normal'% User clicked left-mouse button
                Tempvar.slice_n(4) = round(xy(1,1));
                set(h_taxis,'Value',Tempvar.slice_n(4));
                %set(source,'Selected','on')
            case {'extend','open','alt'} %User did a shift-click
                trial=round(xy(1,1));
                trialnum = sandbox.intrial(trial,3);
                temp = find(sandbox.intrial(:,3)==trialnum);

                if sandbox.intrial(temp,1)==3;
                    sandbox.intrial(temp,1)=1;
                elseif sandbox.intrial(temp,1)==1;
                    sandbox.intrial(temp,1)=3;
                end

                temp = find(sandbox.intrial(:,1)==3);
                idat_red = idat(:,:,:,temp);
                Tempvar.feedb=(['selected trial :' num2str(trialnum)]);
                disp(['selected trial ' num2str(trialnum)]);
                refresh_feedb;
                Tempvar.imagechange = 0;

        end

        refresh_data;
        refresh_image;
        refresh_text;
        Tempvar.imagechange = 0;

    end

    function Click_x_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        xy = get(gca,'CurrentPoint');
        Tempvar.slice_n(2) = round(xy(1,1));
        Tempvar.slice_n(3) = round(xy(1,2));
        Tempvar.slice_l(1:3) = round(Tempvar.slice_n(1:3)./Tempvar.dimratio);
        Tempvar.slice_h(1:3) = round(Tempvar.slice_n(1:3).*Tempvar.dimratio);
        set(h_yaxis,'Value',Tempvar.slice_n(2));
        set(h_zaxis,'Value',Tempvar.slice_n(3));
        refresh_all;
        refresh_text;
        Tempvar.imagechange = 0;
    end

    function Click_y_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        xy = get(gca,'CurrentPoint');
        Tempvar.slice_n(1) = round(xy(1,1));
        Tempvar.slice_n(3) = round(xy(1,2));
        Tempvar.slice_l(1:3) = round(Tempvar.slice_n(1:3)./Tempvar.dimratio);
        Tempvar.slice_h(1:3) = round(Tempvar.slice_n(1:3).*Tempvar.dimratio);
        set(h_xaxis,'Value',Tempvar.slice_n(1));
        set(h_zaxis,'Value',Tempvar.slice_n(3));
        refresh_all;
        refresh_text;
        Tempvar.imagechange = 0;
    end

    function Click_z_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        xy = get(gca,'CurrentPoint');
        Tempvar.slice_n(1) = round(xy(1,2));
        Tempvar.slice_n(2) = round(xy(1,1));
        Tempvar.slice_l(1:3) = round(Tempvar.slice_n(1:3)./Tempvar.dimratio);
        Tempvar.slice_h(1:3) = round(Tempvar.slice_n(1:3).*Tempvar.dimratio);
        set(h_xaxis,'Value',Tempvar.slice_n(1));
        set(h_yaxis,'Value',Tempvar.slice_n(2));
        refresh_all;
        refresh_text;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Image slider callbacks
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Slider_x_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        slider_value_x = get(h_xaxis,'Value');
        Tempvar.slice_n(1) = round(slider_value_x);
        refresh_all;
        Tempvar.imagechange =0;
    end

    function Slider_y_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        slider_value_y = get(h_yaxis,'Value');
        Tempvar.slice_n(2) = round(slider_value_y);
        refresh_all;
        Tempvar.imagechange =0;
    end

    function Slider_z_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        slider_value_z = get(h_zaxis,'Value');
        Tempvar.slice_n(3) = round(slider_value_z);
        refresh_all;
        Tempvar.imagechange =0;
    end

    function Slider_t_Callback(source,eventdata)
        slider_value_t = get(h_taxis,'Value');
        Tempvar.slice_n(4) = round(slider_value_t);
        refresh_data;
        %refresh_figure;
        refresh_image;
        refresh_text;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run SPM5 analysis. Requires SPM5
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function SPM5_Button_Callback(source,eventdata)
        
        fsb_SPM(idat,sandbox);
        idat_spare = idat;
        save idat_temp.mat idat sandbox;
        refresh_map;
        cla;
        refresh_volinfo;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run ICA analysis.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function ICA_Callback(source,eventdata)

        if testmode == 1
            try
                sandbox = fsb_ICA(idat,sandbox);
            catch
                disp('Memory error, running ICA over slices...')
                for x = 1:dim(3)
                    [sandbox icatemp(:,:,x,:)] = fsb_ICA(idat(:,:,x,:),sandbox);
                end
                sandbox.icasig = icatemp;
            end
        else
            str = get(source, 'String');
            val = get(source,'Value');
            sandbox =  fsb_ICA(idat,sandbox,str,val);
        end
        %fsb(sandbox.icasig);
        dummy = sandbox.icasig;
        Map.rc=99;
        Slice_Display_Callback;
        idat_spare = idat;
        save idat_temp.mat idat sandbox;
        refresh_map;
        cla;
        refresh_volinfo;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run fastICA toolbox
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function ICA_Toolbox_Callback(source,eventdata)

        fsb_ICAg(idat,sandbox);
      
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Do 2:1 conversion to increase SNR
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function compress_slices_Callback(source,eventdata)
        dim = size(idat);

        if testmode == 1

            disp(['Old size: ' num2str(dim(1)) ' * '...
                num2str(dim(2)) ' * ' num2str(dim(3))...
                ' Voxels'])

            idat2= idat(:,:,1:2:end,:);
            idat3 = idat(:,:,2:2:end,:);
            idat = idat2+idat3./2;
            dim2 = size(idat);
            disp('Slice compression done')
            disp(['New size: ' num2str(dim2(1)) ' * '...
                num2str(dim2(2)) ' * ' num2str(dim2(3))...
                ' Voxels'])

        else
            return
        end

        idat_spare = idat;
        imean = fsb_calculate_mean(idat);
        save idat_temp.mat idat sandbox;
        refresh_map;
        cla;
        refresh_volinfo;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Keep modified files: Important function!!!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Keep_Modifications
        idat_spare = idat;

        if size(dim,2)>3 && dim(4)>1
            hemodynamics_spare = sandbox.hemodynamics;
            intrial_spare = sandbox.intrial;
            onset_spare = sandbox.onset;
        end;
        refresh_volinfo;
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Load new data
%
% Looks for either .mat or .hdr files and accordingly loads variables
% if .mat is given, it may contain the following fields:
% hemodynamics: A vector of onsets of scans
% onset: A vector containing possible onset responses
% idat: A 3-D or 4-D array containing image data
% hdr: A header for the image data
% intrial: A vector of scans that are part of a trial
% sandbox.discard: A vector of scans that have been deleted
% SPM.mat: A usual SPM.mat file
% if .hdr is given, it will load the image and header provided
% this function requires the use of an external function
% uipickfiles.m which is copyrighted by
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function File_Load_Callback(source,eventdata) %maximum number of files to be loaded at the same time: 63, due to matlab bug

        try

           %[idat sandbox scaninfo] = fsb_load_files([],0,1);
           [idat sandbox scaninfo] = fsb_load_files;
            

            Tempvar.imagechange = 1;

            try
                TR = sandbox.TR;
            catch
                errmsg = lasterror;
                fsb_error_handler(errmsg);
                disp('TR not found, assuming 1 s')
                TR = 1;
                sandbox.TR = TR;
            end
            try
                set(h_xhdr,'String',num2str(sandbox.hdr.dime.pixdim(2)));
                set(h_yhdr,'String',num2str(sandbox.hdr.dime.pixdim(3)));
                set(h_zhdr,'String',num2str(sandbox.hdr.dime.pixdim(4)));
                set(h_xorig,'String',num2str(sandbox.hdr.hist.originator(1)));
                set(h_yorig,'String',num2str(sandbox.hdr.hist.originator(2)));
                set(h_zorig,'String',num2str(sandbox.hdr.hist.originator(3)));
            catch
                errmsg = lasterror;
                fsb_error_handler(errmsg)
                set(h_xhdr,'String',num2str(1));
                set(h_yhdr,'String',num2str(1));
                set(h_zhdr,'String',num2str(1));
            end

            Tempvar.overlayt=0;
            set(gcf,'CurrentAxes',he); cla; axis off;
            set(sandbox_GUI,'Name',['fMRI Sandbox:   ' sandbox.name])

            if isfield(sandbox,'hemodynamics') && ~isempty(sandbox.hemodynamics)
                disp_onset = 1;
            end

            reload_params;
            refresh_volinfo;
            refresh_figure;
            refresh_map;
            refresh_all;
            colormap gray;

        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('ERROR loading files!')
        end
        param.hrf_pred =1;
        set(predsel,'String',num2str(param.hrf_pred));
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Load a high-resolution image for overlay
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Overlay_load_Callback(source,eventdata)

        file_names = uipickfiles('Output','char','REF','.mat|.hdr');
        filnum = size(file_names,1);

        if file_names;
            [pathstr, name, ext, versn] = fileparts(file_names(1,:));

            if filnum==1

                if strcmp(strcat(ext),'.hdr') == 1 || strcmp(strcat(ext),'.img')==1;
                    nidat = avw_img_read_mt(strcat(file_names(1,:)));
                    hidat = int16(nidat.img);
                    clear nidat;
                elseif strcmp(strcat(ext),'.mat') == 1
                    tempdata =load(file_names(1,:));
                    hidat = tempdata.idat;
                end
            end
        end

        dims = size(hidat);

        if dims(3) == 39;
            hidat(:,:,39) = [];
        end
        if length(dims)>3;
            hidat = hidat(:,:,:,1);
        end

        dims = size(hidat);

        if dims(3)>dim(3)
            sc = dims(3)/dim(3);
            if sc ==2;
                hidat(:,:,1:2:end)=[];
            else
                disp('Overlay dimensions do not fit image.')
                disp(['Image slices: ' num2str(dim(3))]);
                disp(['Overlay slices: ' num2str(dims(3))]);

                return
            end
        end

        hidim = size(hidat);
        lodim = size(idat);
        lodim = lodim(1:3);
        Tempvar.dimratio = round(hidim./lodim);
        Tempvar.slice_h=Tempvar.slice_n(1:3).*Tempvar.dimratio;
        Tempvar.loverlay = 1;
        hidat_spare = hidat;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Use a high resolution image to overlay the maps on
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Overlay_use_Callback(source,eventdata)

        if Tempvar.loverlay ~=1
            Overlay_load_Callback
        end

        if (get(h_highres,'Value') == get(h_highres,'Max')) && Tempvar.loverlay ==1;
            Tempvar.hoverlay=1;
            Tempvar.slice_n(1:3) = Tempvar.slice_h(1:3);
        else
            Tempvar.hoverlay =0;
            Tempvar.slice_n = Tempvar.slice_l;
        end

        refresh_data;
        refresh_image;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Change datatype
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function datatype_Callback(source,eventdata)
        str = get(source, 'String');
        val = get(source,'Value');
        what = str{val};

        switch what
            case 'int8'
                idat = int8(idat);
            case 'int16'
                idat = int16(idat);
            case 'single'
                idat = single(idat);
            case 'double'
                idat = double(idat);
        end
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Show spatial transforms
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function transform_display_Callback(source,eventdata)

        if (get(h_transform_display,'Value') == get(h_transform_display,'Max')) ;

            %  Example 1 (3D rotation):
            %	load mri.mat;   old_img = double(squeeze(D));
            for x = 1:dim(4)
                P{x}.mat(1,:) = sandbox.hdr.hist.srow_x;
                P{x}.mat(2,:) = sandbox.hdr.hist.srow_y;
                P{x}.mat(3,:) = sandbox.hdr.hist.srow_z;
                P(x).mat(4,:) = [0 0 0 1];
            end
            
            old_M = P{x}.mat;
            %	old_M = [0.88 0.5 3 -90; -0.5 0.88 3 -126; 0 0 2 -72; 0 0 0 1];
            	new_img = affine(idat, old_M, 2);

            
            
            %             d=[4 0;4 0;4 0];
            %             %function [y1,y2,y3]=coords(p,M1,M2,x1,x2,x3)
            %             % Rigid body transformation of a set of coordinates.
            %             M  = (inv(M2)*inv(spm_matrix(p(1:6)))*M1);
            %             y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*x3 + M(1,4);
            %             y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*x3 + M(2,4);
            %             y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*x3 + M(3,4);
            %
            %
            %             tiny = 5e-2; % From spm_vol_utils.c
            %
            %             Mask = logical(ones(size(y1)));
            %             if ~wrp(1), Mask = Mask & (y1 >= (1-tiny) & y1 <= (dim(1)+tiny)); end;
            %             if ~wrp(2), Mask = Mask & (y2 >= (1-tiny) & y2 <= (dim(2)+tiny)); end;
            %             if ~wrp(3), Mask = Mask & (y3 >= (1-tiny) & y3 <= (dim(3)+tiny)); end;
            %idat = fsb_reslice(idat,sandbox);
        end


        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Save data in memory to disk either as 4-D mat file with hdr field
% or as 3-D analyze files
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function File_Save_Callback(source,eventdata)

        [FileName, PathName] =    uiputfile( ...
            {'*.mat;*.hdr;*.nii','Image Files (*.mat,*.hdr,*.nii)';
            '*.mat','4D Image files (*.mat)'; ...
            '*.hdr','Analyze files (*.hdr)'; ...
            '*.nii','Nifti files (*.nii)'; ...
            '*.*',  'All Files (*.*)'}, ...
            'Save file name');

        if FileName;
            disp(['Saving file ' FileName]);

            [pathstr, name, ext, versn] = fileparts(FileName);

            if ~isempty(ext)
                switch ext
                    case '.mat';
                        fsb_file_write(FileName,idat,sandbox,scaninfo)


                    case '.hdr';

                        idat=full(idat);

                        if isfield(sandbox,'hdr')&& ~isempty(sandbox.hdr);

                            fsb_img_write(idat,name,sandbox.hdr,4,0,scaninfo);
                            sandbox.name = name;
                        else
                            fsb_img_write(idat,name,'template.hdr');
                        end
                    case '.nii'
                        fsb_nii_write(idat,name,sandbox.hdr)
                end
            end
            disp('...done');
        end
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Export data to other formats
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function File_Export_Callback(source,eventdata)

        fsb_export(idat,sandbox)

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This part serves to put in all the functions you might want
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Batch_processing_Callback(source, eventdata)

        %Tempvar.slice_n(1) = 71; Tempvar.slice_n(2) = 92; Tempvar.slice_n(3) = 5; % Center on voxel of interest
        % tic;
        Trial_Sel_All_Callback; % Select all trials
        Trial_Keep_Sel_Callback; % Keep all selected trials
        idat = fsb_realign_3D_ctm(idat,1); % do the realignment with a center of mass algorithm
        Filter_Idat_Spatial_Callback; % smooth volumes
        idat = fsb_normalize(idat,sandbox,3);
        Map_calculate_volume_Callback; % Calculate whole map
        param.hrf_pred = 11; % Select hemodynamics predictor
        Map.corparam = 0.3; % Predefine correlation parameter
        set(h_calcorrelation2_slid,'Value',Map.corparam); % Set threshold slider to correct position
        set(predsel, 'Value',param.hrf_pred); % Set correct predictor
        Tempvar.mapdisp = 1; % Map output on
        set(h_mapdisp,'Value',1); %Set Map display radiobutton to pressed
        refresh_all;
        %t= toc;
        %disp(['Batch processing image data took ' num2str(t) ' seconds'])

    end

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Get Design information
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Design_load_Callback(source, eventdata)

        str = get(source, 'String');
        val = get(source,'Value');
        what = str{val};
        dim = size(idat);

        sandbox = fsb_get_hemodynamics(dim,sandbox,what);

        Tempvar.imagechange = 1;

        [sandbox prms] = fsb_import_sandbox(sandbox,idat);


        if size(dim,2)>3 && dim(4) ~= size(sandbox.hemodynamics,1) || max(idat(:))==1;
            dim(4) = size(sandbox.hemodynamics,1);

            bhemodynamics = (sandbox.hemodynamics+0.2)*(dim(2)/4);
            idat = ones(dim(1),dim(2),dim(3),dim(4));
            idat_spare = idat;

        end
        disp_onset =1;
        param.hrf_pred =1;
        %refresh_volinfo;
        reload_params;
        refresh_map;
        refresh_all;
        try
            disp(['Hemodynamics: '  num2str(size(sandbox.hemodynamics,1))]);
            disp(['Intrial: '  num2str(size(sandbox.intrial,1))]);
            disp(['Onset: '  num2str(size(sandbox.onset,1))]);
            disp(['Size (d) ': num2str(dim(4))]);
            disp(['Successful trials: ' num2str(sandbox.origtrial)]);
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
        end

        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Save mean image to disk
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Write_Mean_Callback(source,eventdata)
        %[file_name,path_name] = uiputfile('subject','Save file name');
        [file_name, PathName] =    uiputfile( ...
            {'*.mat;*.hdr;*.nii','Image Files (*.mat,*.hdr,*.nii)';
            '*.mat','4D Image files (*.mat)'; ...
            '*.hdr','Analyze files (*.hdr)'; ...
            '*.nii','Nifti files (*.nii)'; ...
            '*.*',  'All Files (*.*)'}, ...
            'Save file name');
        [pathstr, name, ext, versn] = fileparts(file_name);

        if ~isempty(ext)
            switch ext
                case '.mat';
                    fsb_file_write(file_name,imean,sandbox,scaninfo)


                case '.hdr';

                    idat=full(idat);

                    if isfield(sandbox,'hdr')&& ~isempty(sandbox.hdr);

                        fsb_img_write(imean,name,sandbox.hdr,4,0,scaninfo);
                        sandbox.name = name;
                    else
                        fsb_img_write(imean,name,'template.hdr');
                    end
                case '.nii'
                    fsb_nii_write(imean,name,sandbox.hdr)
            end
        end
        %         if isfield(sandbox,'hdr')& size(sandbox.hdr)>0;
        %             fsb_img_write(imean,file_name_analyze,sandbox.hdr,scaninfo,0,4);
        %         else
        %             fsb_img_write(imean,file_name_analyze,'template.hdr');
        %         end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Save selected image to disk
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Write_Selected_Callback(source,eventdata)
        [file_name_analyze,path_name] = uiputfile('subject','Save file name');

        if isfield(sandbox,'hdr')& size(sandbox.hdr)>0;
            fsb_img_write(idat(:,:,:,Tempvar.slice_n(4)),file_name_analyze,sandbox.hdr,scaninfo,0,4);
        else
            fsb_img_write(idat(:,:,:,Tempvar.slice_n(4)),file_name_analyze,'template.hdr');
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Threshold with snr map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Threshold_SNR_Slider_Callback(source, eventdata)

        Tempvar.threshold_SNR = get(source,'Value');
        idat = idat_spare;
        dim = size(idat);
        noise = single(idat(1:2,1:round(dim(2)/2),round(dim(3)/2):end,:));
        noise = std(single(noise),0,4);
        noise_std = mean(noise(:));
        intensity_threshold = noise_std*Tempvar.threshold_SNR;
        disp(['Intensity threshold: ' num2str(intensity_threshold)]);

        idat = fsb_denoise_idat(idat_spare,intensity_threshold,0);
        disp(['SNR thresholding done,threshold at SNR = ' num2str(Tempvar.threshold_SNR)]);
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Write cluster timecourse
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function TC_write_Cluster_Callback(source,eventdata)

        figure(114);clf;
        otto = plot(trial_hrf2_normalized, 'd','MarkerSize',5);
        hold on;
        ottocal = hsv(size(otto,1));

        for x=1:size(otto)
            set(otto(x),'MarkerEdgeColor',ottocal(x,:));
            set(otto(x),'MarkerFaceColor',ottocal(x,:));
        end

        plot(trial_normalized,'k','LineWidth',3);
        plot(hem2trial_normalized, ':*k','LineWidth',2);

        psc_string = sprintf('%1.2f',Tempvar.psc);
        text(0.02,0.08,['signal change = ' psc_string ' %'],...
            'Interpreter','none','FontSize',8,'Units','normalized','Color','Black','BackgroundColor','White');

        xlabel('trial timepoint');
        ylabel('% signal change');
        Title('Signal change over trials');
        %legend('Timepoints','Mean trial change','Predicted trial change')

        [FileName, PathName] =    uiputfile( ...
            {'*.mat','Image Files (*.mat)'},...
            'Pick a file');

        if FileName
            disp(['Saving file ' FileName]);
            TC_cluster = trial_hrf2_normalized;
            save (FileName,'TC_cluster')
        end
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Select ROI size
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function ROI_size_Slider_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        Tempvar.roi = floor(get(h_roi_slid,'Value'));
        Tempvar.feedb =['ROI size : ' num2str((2*Tempvar.roi)^2) ' voxel'];

        refresh_all;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Delete ROI contents - useful if you want to remove tissue bridges before
% extracting brains
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function ROI_voxels_delete_Callback(source,eventdata)
        idat(Tempvar.roiind(1,1):Tempvar.roiind(1,2),...
            Tempvar.roiind(2,1):Tempvar.roiind(2,2),Tempvar.roiind(3,1):Tempvar.roiind(3,2),:) = 0;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Keep ROI contents - useful if you want to keep just parts of the volume
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function ROI_voxels_keep_Callback(source,eventdata)
        idat = idat(Tempvar.roiind(1,1):Tempvar.roiind(1,2),...
            Tempvar.roiind(2,1):Tempvar.roiind(2,2),Tempvar.roiind(3,1):Tempvar.roiind(3,2),:);
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Extract brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Extract_Slider_Callback(source,eventdata)
        Tempvar.slice_n(4) =1;
        Tempvar.extr_param = floor(get(h_extrbrain_slid,'Value'));
        b_idat = fsb_get_extraction_mask(imean,imean_avg,Tempvar.extr_param);
        idat = b_idat;
        Tempvar.feedb =['Extraction threshold ' num2str(Tempvar.extr_param) ];
        refresh_figure;
        refresh_data;
        refresh_image;
        idat = idat_spare;
    end

    function Brain_Extract_Button_Callback(source,eventdata)
        pack;
        Tempvar.slice_n(4) =1;
        h = waitbar(0,'Extracting brain...');
        [dim(1),dim(2),dim(3),dim(4)] = size(idat_spare);
        waitbar(1/10);
        ithres = b_idat==0;
        ithres2 = int32(find(b_idat==0));
        idat=reshape(idat_spare,dim(1)*dim(2)*dim(3),dim(4));%reshape to be able to use linear index
        waitbar(2/10);
        idat(ithres2,:)=NaN; %kick out voxels of no interest
        idat=reshape(idat,dim(1),dim(2),dim(3),dim(4)); %put everything back in shape

        if size(dim,2)>3 && dim(4)>1;
            ithres_map = repmat(ithres,[1 1 1 size(sandbox.hemodynamics,2)]);
            sandbox.stats.r_cor(ithres_map) = NaN;
            sandbox.stats.brob(ithres_map) = NaN;
            sandbox.stats.b_est(ithres_map) = NaN;
            sandbox.stats.p_rob(ithres_map) = NaN;
            sandbox.stats.rob_est(ithres_map) = NaN;
            sandbox.stats.t_test(ithres_map) = NaN;
            r_cor_red = zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
            r_cor_neg =zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
        end

        waitbar(6/10)
        Tempvar.feedb ='Extraction finished' ;
        idat_spare = idat;
        waitbar(7/10)
        refresh_volinfo;
        waitbar(10/10)
        close(h);
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Detect artefacts
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Artefact_detection_Slider_Callback(source,eventdata)%#ok<INUSD> % detect artefacts

        Tempvar.detart_param = get(h_detart_slider,'Value');
        sandbox.discard =[];
        sandbox.discard  = fsb_det_artifact_mod(idat,Tempvar.detart_param);

        if size(sandbox.hemodynamics,2)>20; % Change this to 9 if you want to include motion sensor data
            motionstd1 = std(sandbox.hemodynamics(:,8));
            motionstd2 = std(sandbox.hemodynamics(:,9));
            motionmean1 = mean(sandbox.hemodynamics(:,8));
            motionmean2 = mean(sandbox.hemodynamics(:,9));
            discardjaw = find(sandbox.hemodynamics(:,8)>(motionmean1+motionstd1*1.5));
            discardpo = find(sandbox.hemodynamics(:,9)>(motionmean2+motionstd2*1.5));
            discardmotion = union (discardjaw, discardpo);
            sandbox.discard = union(sandbox.discard,discardmotion);
        end

        [sandbox Tempvar] = display_discard(sandbox,Tempvar);

        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Remove artefacts
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Artefact_Remove_Button_Callback(source,eventdata) % remove artefacts
        if size(sandbox.discard,1) == dim(4)
            disp('Trying to remove entire dataset, aborting')
            Tempvar.feedb = 'Trying to remove entire dataset, aborting';
            refresh_feedb;
            return;
        end

        tic;
        Tempvar.imagechange = 1;
        h = waitbar(0, 'Removing artefacts...');
        waitbar(1/10);
        hemodynamics_spare = sandbox.hemodynamics;
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        sandbox.artefacts = sandbox.hemodynamics;

        sandbox = fsb_recalc_onsets(sandbox,sandbox.discard);

        for x = 1:size(sandbox.discard,2);
            sandbox.artefacts(x,size(sandbox.hemodynamics+x))=1;
        end

        if size(sandbox.hemodynamics,1)== dim(4);
            sandbox.hemodynamics(sandbox.discard,:)=[];
        end
        if size(sandbox.onset,1)== dim(4);
            sandbox.onset(sandbox.discard,:)=[];
        end
        waitbar(2/10)
        if size(sandbox.intrial,1)== dim(4);
            sandbox.intrial(sandbox.discard,:)=[];
            Tempvar.intrial_nr =[];
            for x = 1:size(sandbox.intrial,2);
                Tempvar.intrial_nr(:,x) = bwlabeln(sandbox.intrial(:,x));
            end
        end
        waitbar(4/10)
        idat(:,:,:,sandbox.discard) = [];
        sandbox.discard = zeros(dim(4),1);
        Tempvar.disp_disc = 0;
        Tempvar.do=dim(4);
        dim = size(idat);
        if Tempvar.slice_n(4) > dim(4);
            Tempvar.slice_n(4) = 1;
        end
        waitbar(5/10)
        t = toc;
        Tempvar.feedb =['Artifacts removed in ' num2str(t) ' seconds'];
        trialmax = max(Tempvar.intrial_nr(:,3));

        for x = 1: trialmax;
            test = find(Tempvar.intrial_nr(:,3) == x);
            trial_length(x) = size(test,1);
        end
        try
            disp (['Trial length: ' num2str(trial_length)]);
            disp (['Number of trials: ' num2str(trialmax)]);
        catch
            disp('All trials removed');
        end

        sandbox.discard = [];

        idat_spare =idat;
        hemodynamics_spare = sandbox.hemodynamics;
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        waitbar(8/10)
        close(h);
        refresh_map;
        refresh_volinfo;
        refresh_all;
        Tempvar.imagechange = 0;
        sandbox.fullscan = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% remove selected defective volume
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Delete_Selected_Callback(source,eventdata)
        tic;
        Tempvar.imagechange = 1;
        hemodynamics_spare = sandbox.hemodynamics;
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        discardv = Tempvar.slice_n(4);
        sandbox = fsb_recalc_onsets(sandbox,discardv);
        if size(sandbox.hemodynamics,1)== dim(4);
            sandbox.hemodynamics(discardv,:)=[];
        end
        if size(sandbox.discard,1)>1
            sandbox.discard(discardv==sandbox.discard,:)=[];
        end
        if size(sandbox.onset,1)== dim(4);
            sandbox.onset(discardv,:)=[];
        end
        try
        if size(sandbox.sensdata,1)== dim(4);            
            sandbox.sensdata(discardv,:)=[]; 
        end
        catch
            disp('No sensor data field present')
        end
        if size(sandbox.intrial,1)== dim(4);
            sandbox.intrial(discardv,:)=[];
            Tempvar.intrial_nr =[];
            for x = 1:size(sandbox.intrial,2);
                Tempvar.intrial_nr(:,x) = bwlabeln(sandbox.intrial(:,x));
            end
        else
            Tempvar.intrial_nr = [];
        end
        idat(:,:,:,discardv) = [];
        %Tempvar.disp_disc = 0;
        Tempvar.do=dim(4);
        dim = size(idat);
        if size(Tempvar.slice_n,1)<4 || Tempvar.slice_n(4) > dim(4);
            Tempvar.slice_n(4) = 1;
        end

        t = toc;
        Tempvar.feedb =['Artifact removed in ' num2str(t) ' seconds'];

        try
            trialmax = max(Tempvar.intrial_nr(:,3));
            trial_length = [];
            for x = 1: trialmax;
                test = find(Tempvar.intrial_nr(:,3) == x);
                trial_length(x) = size(test,1);
            end
            disp (['Trial length: ' num2str(trial_length)]);
            disp (['Number of trials: ' num2str(trialmax)]);
        catch
            disp('No trials defined yet')
        end

        if isfield(sandbox,'done.volumedelete')
            sandbox.done.volumedelete = [sandbox.done.volumedelete discardv];
        else
            sandbox.done.volumedelete = discardv;
        end
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        hemodynamics_spare = sandbox.hemodynamics;
        refresh_feedb;
        refresh_volinfo;
        refresh_map;
        refresh_all;
        Tempvar.imagechange = 0;
        sandbox.fullscan = 0;
        Keep_Modifications;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% add one baseline volume in intrial file before every trial to extend the baseline
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_add_baseline_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        tic;

        intrial2 = sandbox.intrial;
        sandbox.intrial(:,3) = bwlabeln(intrial2(:,1)); % Label each trial with own number

        if max(sandbox.intrial(:,3))>1 && sandbox.intrial(1,1)==0;

            for x = 1:max(sandbox.intrial(:,3));
                scanind = find (sandbox.intrial(:,3)==x);
                scanind = [min(scanind)-1; scanind];
                sandbox.intrial(scanind,1)=1;
            end;
            Tempvar.feedb = ['1 pre-Baseline Volume added to trial' ];
            t = toc;
        else
            Tempvar.feedb =['Only one trial loaded or no space left'];
        end
        intrial2 = sandbox.intrial;
        sandbox.intrial(:,3) = bwlabeln(intrial2(:,1)); % Label each trial again

        if isfield(sandbox,'done.baselineadd')
            sandbox.done.baselineadd = sandbox.done.baselineadd+1;
        else
            sandbox.done.baselineadd = 1;
        end

        refresh_feedb;
        refresh_map;
        refresh_all;
        Tempvar.imagechange =0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% remove one baseline volume in intrial file before every trial
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_rem_baseline_Callback(source,eventdata)
        Tempvar.imagechange = 1;

        intrial2 = sandbox.intrial;
        intrial2(intrial2(:,2)==1,3)=0;
        sandbox.intrial(:,3) = bwlabeln(intrial2(:,3)); % Label each trial with own number

        if max(sandbox.intrial(:,3))>1 && sandbox.intrial(1,1)==0;

            for x = 1:max(sandbox.intrial(:,3));
                scanind = find (sandbox.intrial(:,3)==x);
                sandbox.intrial(scanind(1),1)=0;
                sandbox.intrial(scanind(1),3)=0;
            end;
            Tempvar.feedb = ['1 pre-Baseline Volume removed from trial' ];
        else
            Tempvar.feedb =['Only one trial loaded or no space left'];
            sandbox.intrial(sandbox.intrial(:,3)==0,1)=0;
        end

        if isfield(sandbox,'done.baselineremove')
            sandbox.done.baselineremove = sandbox.done.baselineremove+1;
        else
            sandbox.done.baselineremove = 1;
        end

        intrial2 = sandbox.intrial;
        sandbox.intrial(:,3) = bwlabeln(intrial2(:,1)); % Label each trial again
        refresh_feedb;
        refresh_map;
        refresh_all;
        Tempvar.imagechange =0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Select/Unselect all trials
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_Sel_All_Callback(source, eventdata)

        sandbox.intrial(sandbox.intrial(:,1)==1)=3;
        trial_nr = max(sandbox.intrial(:,3));
        disp(['Number of selected trials: ' num2str(trial_nr)]);
        refresh_image;
        refresh_feedb;

    end

    function Trial_Sel_None_Button_Callback(source, eventdata)

        sandbox.intrial(sandbox.intrial(:,1)==3)=1;
        refresh_image;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Remove intrial and hemodynamics indexes for selected trials from
% sandbox.intrial/hemodynamics file
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_label_selected_Callback(source, eventdata)

        sandbox.hemodynamics(sandbox.intrial(:,1)==3,:)=0;
        sandbox.intrial(sandbox.intrial(:,1)==3,1:3)=0;
        refresh_map;
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Put an artefact label on the selected volume
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_label_selected_Callback(source, eventdata)

        sandbox.discard = [sandbox.discard; Tempvar.slice_n(4)];
        [sandbox Tempvar] = display_discard(sandbox,Tempvar);
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% keep only selected  trials
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_Keep_Sel_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        tic;

        intrial_select = find(sandbox.intrial(:,1) ==3);
        discard = find(sandbox.intrial(:,1)<3);

        sandbox = fsb_recalc_onsets(sandbox,discard);

        if isempty(intrial_select);
            disp('No trials selected');
            return;
        end

        sandbox.done.trialselect = intrial_select;


        if size(sandbox.hemodynamics,1)== dim(4);
            sandbox.hemodynamics=sandbox.hemodynamics(intrial_select,:);
        end

        if size(sandbox.onset,1)== dim(4);
            sandbox.onset=sandbox.onset(intrial_select,:);
        end

        if size(sandbox.sensdata,1)== dim(4);
            sandbox.sensdata=sandbox.sensdata(intrial_select,:);
        end

        if size(sandbox.intrial,1)== dim(4);
            sandbox.intrial=sandbox.intrial(intrial_select,:);
            Tempvar.intrial_nr =[];

            for x = 1:size(sandbox.intrial,2);
                Tempvar.intrial_nr(:,x) = bwlabeln(sandbox.intrial(:,x));
            end

        end

        idat = idat(:,:,:,intrial_select);
        Tempvar.do=dim(4);
        dim = size(idat);

        if Tempvar.slice_n(4) > dim(4);
            Tempvar.slice_n(4) = 1;
        end

        r_cor_spare =sandbox.stats.r_cor;
        r_cor_red_spare = r_cor_red;
        r_cor_neg_spare =r_cor_neg;
        hemodynamics_spare = sandbox.hemodynamics;

        intrial = sandbox.intrial;
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        idat_spare = idat;
        sandbox.intrial(:,1) =1;
        sandbox.fullscan = 0;

        t = toc;
        disp(['Selected volumes cropped in ' num2str(t) ' seconds']);
        refresh_volinfo;
        refresh_map;
        refresh_all;
        Tempvar.imagechange = 0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% delete selected  trials
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_Del_Sel_Callback(source,eventdata)

        intrial_select = find(sandbox.intrial(:,1) ==3);

        if size(intrial_select,1) == dim(4)
            disp('Trying to remove entire dataset, aborting')
            Tempvar.feedb = 'Trying to remove entire dataset, aborting';
            refresh_feedb;
            return;
        end
        if isfield(sandbox,'done.trialremove')
            sandbox.done.trialremove = [sandbox.done.trialremove intrial_select];
        else
            sandbox.done.trialremove = intrial_select;
        end
        Tempvar.imagechange = 1;
        tic;
        r_cor_spare =sandbox.stats.r_cor;
        r_cor_red_spare = r_cor_red;
        r_cor_neg_spare =r_cor_neg;

        hemodynamics_spare = sandbox.hemodynamics;
        intrial_spare = sandbox.intrial;
        onset_spare = sandbox.onset;
        idat_spare = idat;


        sandbox.discard = union(intrial_select,sandbox.discard);
        sandbox.discard = unique(sandbox.discard);

        sandbox = fsb_recalc_onsets(sandbox,sandbox.discard);

        try
            sandbox.hemodynamics(intrial_select,:) =[];
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('No hemodynamics file loaded');
        end

        try
            sandbox.sensdata(intrial_select,:) = [];
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('No sensor information present');
        end

        try
            hemodynamics_deriv(intrial_select,:) =[];
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('No hemodynamics_deriv file loaded');
        end

        try
            sandbox.discard(intersect(sandbox.discard,intrial_select),:) =[];
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            sandbox.discard = zeros(dim(4),1);
            disp('Discard file not concatenated');
        end

        try
            ctr_mass(intrial_select,:) =[];
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('No ctr_mass file loaded');
        end

        if size(sandbox.onset,1)== dim(4);
            sandbox.onset(intrial_select,:) =[];
        end

        if size(sandbox.intrial,1)== dim(4);
            sandbox.intrial(intrial_select,:) =[];
            Tempvar.intrial_nr =[];

            for x = 1:size(sandbox.intrial,2);
                Tempvar.intrial_nr(:,x) = bwlabeln(sandbox.intrial(:,x));
            end

        end

        idat(:,:,:,intrial_select) =[];
        Tempvar.do=dim(4);
        dim = size(idat);
        if Tempvar.slice_n(4) > dim(4);
            Tempvar.slice_n(4) = 1;
        end
        t = toc;
        Tempvar.feedb =['Selected volumes deleted in ' num2str(t) ' seconds'];

        trialmax = max(Tempvar.intrial_nr(:,3));
        trial_length = zeros(1,trialmax);

        for x = 1: trialmax;
            test = find(Tempvar.intrial_nr(:,3) == x);
            trial_length(x) = size(test,1);
        end
        if trialmax == 0;
            trial_length = 0;

        end
        disp (['Trial length: ' num2str(trial_length)]);
        disp (['Number of trials: ' num2str(trialmax)]);
        idat_spare = idat;
        intrial_spare = sandbox.intrial;
        sandbox.fullscan = 0;

        refresh_volinfo;
        refresh_map;
        refresh_all;
        Tempvar.imagechange = 0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Remove empty space around the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Crop_Callback(source,eventdata)
        tic;
        idat_spare = idat;
        pack;
        [idat sandbox] = fsb_crop_brain(idat,sandbox);
        idat_spare = idat;
        [Tempvar.slice_n(1) Tempvar.slice_n(2) Tempvar.slice_n(3)] = size(idat);
        t = toc;

        Tempvar.feedb =['Cropping done in ' num2str(t) ' seconds'];

        refresh_volinfo;
        refresh_map;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Pad it with zeros to get predefined image format
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Pad_Callback(source,eventdata)
        tic;
        idat = fsb_pad_idat(idat,0);
        % idat = p_idat;
        idat_spare = idat;
        t = toc;
        Tempvar.feedb =['Padding done in ' num2str(t) ' seconds'];
        cla;
        refresh_volinfo;
        refresh_map;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Do analyses on resting state periods
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Baseline_conversion_Callback(source,eventdata)

        str = get(source, 'String');
        val = get(source,'Value');
        Tempvar.imagechange = 1;

        switch str{val};

            case 'Convert trial structure'

                disp('Converting trial structure to also include aborted trials...');
                sandbox = fsb_get_all_trial(sandbox);

            case 'Interpolate artefacts and trial data'

                disp('Interpolating artefacts and trials run by run...');
                sandbox.filttype = 'Stop';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;
                filtopts.discard = 1;
                filtopts.invert = 1;
                filtopts.filter =1;
                filtopts.mintrial = 3;
                [idat sandbox] = fsb_interpolate_filter(idat,sandbox,filtopts);

            case 'Create regressor for current ROI'

                param.hrf_pred = 90;
                sandbox.hemodynamics(:,90) = Tempvar.roimean;

                set(predsel,'String',num2str(param.hrf_pred));
                disp ('defaulting to ROI predictor');
                Tempvar.feedb = ('defaulting to ROI predictor');

                refresh_map;
                refresh_image;

        end
        refresh_all
        Tempvar.imagechange = 0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Normalize images
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Normalize_Callback(source,eventdata)

        str = get(source, 'String');
        val = get(source,'Value');

        tic;
        Tempvar.imagechange = 1;
        idat_spare = idat;

        switch str{val};
%             case 'Normalization Toolbox' % Not yet implemented
%                 disp('Starting Toolbox');
%                 idat = fsb_normalize_GUI(idat,sandbox,str,val);
            case 'Remove outliers'
                idat = fsb_remove_outliers(idat,imean);
                sandbox.done.outliercorrection = 1;
                disp('Removing outliers');
            case 'Normalize Volumes'
                idat = fsb_normalize_scans(idat,0);
                sandbox.done.volume_normalization = 1;
                disp('Normalizing volumes with average intensity');
            case 'Detrend scans'
                disp('Correcting scans for linear trends');
                idat = fsb_detrend_scan(idat,sandbox,0,0);
                sandbox.done.detrending_linear = 1;
            case 'Detrend scans nonlinear'
                disp('Correcting scans for nonlinear trends');
                idat = fsb_detrend_scan(idat,sandbox,0,1);
                sandbox.done.detrending_nonlinear = 1;
            case 'Interpolate Artefacts'
                disp('Interpolating artefacts parallel run by run...');
                idat = fsb_interpolate_parallel_scans(idat,sandbox);
                sandbox.done.artefact_interpolation = 1;
            case 'Interpolate Artefacts and band-pass filter data'
                disp('Interpolating artefacts run by run and applying filter...');
                sandbox.filttype = 'Stop';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;
                idat = fsb_interpolate_filter_scans(idat,sandbox);
                sandbox.done.interpolate_bandpass = sandbox.filt_freq;
            case 'Interpolate Artefacts and high-pass filter data'
                disp('Interpolating artefacts run by run and applying filter...');
                sandbox.filttype = 'High';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;
                idat = fsb_interpolate_filter_scans(idat,sandbox);
                sandbox.done.interpolate_bandpass = sandbox.filt_freq;
            case 'Interpolate Artefacts and low-pass filter data'
                disp('Interpolating artefacts run by run and applying filter...');
                sandbox.filttype = 'Low';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;
                idat = fsb_interpolate_filter_scans(idat,sandbox);
                sandbox.done.interpolate_bandpass = sandbox.filt_freq;
            case 'Enhance contrast'
                disp('Improving gray-white matter contrast');
                idat = fsb_enhance_contrast(idat,sandbox);
                sandbox.done.contrast_enhancement = 1;
            case 'Normalize Sessions'
                idat = fsb_normalize(idat,sandbox,5);
                sandbox.done.session_normalization = 1;
                disp('Normalizing sessions with average intensity');
            case 'Normalize Runs'
                idat = fsb_normalize(idat,sandbox,4);
                sandbox.done.run_normalization = 1;
                disp('Normalizing runs with average intensity');
            case 'Normalize Trials'
                idat = fsb_normalize(idat,sandbox,3);
                sandbox.done.trial_normalization = 1;
                disp('Normalizing trials with average intensity');
            case 'Normalize Baseline multiplicative'
                idat = fsb_normalize_baseline(idat,sandbox);
                sandbox.done.baseline_normalization_multiplicative = 1;
                disp('Normalizing trials with their baseline intensity - multiplicative');
            case 'Normalize Baseline divisive'
                idat = fsb_normalize_baseline_divisive(idat,sandbox);
                sandbox.done.baseline_normalization_divisive = 1;
                disp('Dividing trials through their baseline intensity');
            case 'Normalize Baseline additive'
                idat = fsb_normalize_baseline_add(idat,sandbox,param.hrf_pred);
                sandbox.done.baseline_normalization_additive = 1;
                disp('Normalizing trials with their baseline intensity - additive');
            case 'Normalize Baseline Z-score'
                idat = fsb_normalize_baseline_zscore(idat,sandbox);
                sandbox.done.baseline_normalization_zscores = 1;
                disp('Normalizing trials with their baseline intensity - z-scores');
            case 'Normalize Z-score'
                idat = fsb_normalize_zscore(idat);
                sandbox.done.volume_normalization_zscores = 1;
                disp('Normalizing volumes to z-scores');
            case 'Normalize Noise'
                idat = fsb_normalize_noise(idat,sandbox);
                sandbox.done.normalization_noise_out_of_brain = 1;
                disp('Normalizing trials with out-of-brain noise');
            case 'Normalize ROI'
                idat = fsb_normalize_ROI(idat,Tempvar.roimean);
                sandbox.done.normalization_ROI = 1;
                sandbox.done.normalization_ROImean = Tempvar.roimean;
                disp('Normalizing trials with ROI informatioin');
        end

        try
            idat(isnan(idat))=0;
            idat(idat<0) = 0;
            idat(idat>32500) = 32500;
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('Memory problems, skipping thresholding.. ');
        end
        imean = fsb_calculate_mean(idat);
        idat_spare = idat;

        refresh_map;
        refresh_volinfo;
        refresh_all;
        Tempvar.imagechange = 0;
        t = toc;
        disp(['Normalization done in ' num2str(t) ' seconds']);

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Header handling
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HDR_Callback(source,eventdata)
        hdr_x = get(h_xhdr,'String');
        hdr_y = get(h_yhdr,'String');
        hdr_z = get(h_zhdr,'String');
        orig_x = get(h_xorig,'String');
        orig_y = get(h_yorig,'String');
        orig_z = get(h_zorig,'String');
        set(h_xhdr,'String',num2str(hdr_x));
        set(h_yhdr,'String',num2str(hdr_y));
        set(h_zhdr,'String',num2str(hdr_z));
        set(h_xorig,'String',num2str(orig_x));
        set(h_yorig,'String',num2str(orig_y));
        set(h_zorig,'String',num2str(orig_z));
        sandbox.hdr.dime.pixdim(2) = str2double(hdr_x);
        sandbox.hdr.dime.pixdim(3) = str2double(hdr_y);
        sandbox.hdr.dime.pixdim(4) = str2double(hdr_z);
        sandbox.hdr.hist.originator(1) = str2double(orig_x);
        sandbox.hdr.hist.originator(2) = str2double(orig_y);
        sandbox.hdr.hist.originator(3) = str2double(orig_z);
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Diagnostics
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Slice_Display_Callback(source, eventdata)

        if Tempvar.mapcalc == 0 && Tempvar.mapdisp == 1;
            Map_calculate_volume_Callback
        end

        try
            str = get(source, 'String');
            val = get(source,'Value');
            inparg.maptype = str{val};
            Map.maptype = inparg.maptype;
        catch
            inparg.maptype = 'Plot slices';
            Map.maptype = inparg.maptype;
        end

        inparg.maptype = Map.maptype;
        inparg.hoverlay = Tempvar.hoverlay;
        inparg.calcorr = Tempvar.mapdisp;
        inparg.dummy = dummy;
        inparg.slicez_data = slicez_data;
        inparg.rc = Map.rc;
        inparg.tsfx = Map.tsfx;
        inparg.tsfy = Map.tsfy;
        inparg.shiftx = Map.shiftx;
        inparg.shifty = Map.shifty;
        inparg.corparam = Map.corparam;
        inparg.probparam = Map.probparam;
        inparg.corrstring = get(h_probability_button,'String');
        inparg.probstring = get(h_probability_slid,'String');
        inparg.hrf_pred = param.hrf_pred;
        inparg.newmap = newmap;
        inparg.slice_n = Tempvar.slice_n;
        inparg.roi = Tempvar.roi;
        inparg.hnoise = 6;
        inparg.hrf_pred = param.hrf_pred;
        inparg.nscal = 0;
        inparg.idat = idat; % this is done just for testing
        inparg.imean = imean; % this is done just for testing

        if Tempvar.hoverlay ==1;
            inparg.hidat = hidat;
        end

        fsb_diag(idat,imean,sandbox,inparg);
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Overlay different trials in subplot
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% Turn trial overlays on

    function Trial_display_timecourse_Button_Callback(source,eventdata)

        Tempvar.imagechange = 1;
        if (get(h_overlaytrial,'Value') == get(h_overlaytrial,'Max')) ;
            Tempvar.overlayt=1;
            Tempvar.feedb ='Trial overlay on ';
            set(gcf,'CurrentAxes',he); cla; axis on;
        else
            Tempvar.overlayt = 0;
            Tempvar.feedb ='Trial overlay off ';
            set(gcf,'CurrentAxes',he); cla; axis off;
        end

        refresh_all;
        Tempvar.imagechange = 0;
    end

% Turn cluster overlays on

    function Trial_TC_voxel_cluster_Callback(source,eventdata)

        Tempvar.imagechange = 1;
        if (get(h_overlayclus,'Value') == get(h_overlayclus,'Max')) ;
            clustselect = 1;
            Tempvar.feedb ='Cluster timecourse display on ' ;

        else
            clustselect = 0;
            Tempvar.feedb ='ROI timecourse display on ' ;

        end
        refresh_all;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate cluster data overlays
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_Overlay_Callback(source,eventdata)

        if Tempvar.imagechange ==1;

            if analtrial ==1;
                hemodynamics_spare = sandbox.hemodynamics;
                intrial_spare = sandbox.intrial;
                sandbox.hemodynamics = hemodynamics_red;
                trial_index = sandbox.intrial(:,1)==3;
                sandbox.intrial = sandbox.intrial(trial_index,:);
                idat_spare = idat;
                idat = idat_red;
            end

            tslice = Tempvar.slice_l(1:3);

            if     clustselect==1 && r_cor_red(tslice(1),tslice(2),tslice(3),param.hrf_pred)>0;

                [B] = FloodFill3D_s(r_cor_red(:,:,:,param.hrf_pred), [tslice(1) tslice(2) tslice(3)]);

                if (sum(sum(sum(B))))>1;
                    cluster_size = sum(sum(sum(B)));
                    B = repmat(B,[1 1 1 size(idat,4)]);
                    clust_hrf= idat(B>0);
                    clustc2=reshape(clust_hrf,[],size(idat,4));
                    trial_hrf = squeeze(nanmean(clustc2));
                    trial_hrf = trial_hrf';
                else
                    trial_hrf= squeeze(idat(tslice(1),tslice(2),tslice(3),:));
                    cluster_size = 1;
                end

            elseif clustselect==0 || clustselect==1 && r_cor_red(tslice(1),tslice(2),tslice(3),param.hrf_pred)==0;

                try
                    trial_hrf = Tempvar.roimean;
                    cluster_size = size(find(Tempvar.roimask>0),1);
                catch
                    errmsg = lasterror;
                    fsb_error_handler(errmsg);
                    disp('ROI edge exceeds matrix dimensions, aborting operation')
                    if analtrial ==1;
                        sandbox.hemodynamics = hemodynamics_spare;
                        sandbox.intrial = intrial_spare;
                        idat = idat_spare;
                    end
                    return;
                end

            end

            % Tempvar.do cluster calculations with subfunction
            tc.ROI_data.Y =  trial_hrf;
            tc.ROI_data.block_rows{1} = find(ones(size(trial_hrf,1),1)>0);
            tc.ROI_data.regions{1}.Y = trial_hrf;
            tc.sandbox = sandbox;
            tc.sandbox.sb{1}.intrial = sandbox.intrial;
            tc.sandbox.sb{1}.hemodynamics = sandbox.hemodynamics;
            tc.sandbox.sb{1}.onset = sandbox.onset;

            tc.hrf_pred = param.hrf_pred;
            tc.display = 1;
            tc.reshuffle = 0;
            try
                tcm = fsb_calc_cluster_tc(tc,0);
            catch
                errmsg = lasterror;
                fsb_error_handler(errmsg);
                disp('Issues with cluster algorithm, aborting...')
                if analtrial ==1;
                    sandbox.hemodynamics = hemodynamics_spare;
                    sandbox.intrial = intrial_spare;
                    idat = idat_spare;
                end
                return;
            end
            set(0,'CurrentFigure',sandbox_GUI);

            trial_hrf2_normalized = tcm.trial_hrf2_normalized{1};
            trial_normalized = tcm.trial_normalized{1};
            Tempvar.psc = tcm.psc{1};
            snr_trial = tcm.snr_trial{1};
            corr_trial = tcm.corr{1};
            hem2trial_normalized = tcm.session_hrf_normalized{1};

            set(sandbox_GUI,'CurrentAxes',he); cla; axis on; hold on;
            cla;

            if ~isfield(sandbox,'plotall')
                sandbox.plotall = 0;
            end

            if sandbox.plotall == 1;
                otto = plot(trial_hrf2_normalized, 'd','MarkerSize',5);
                ottocal = hsv(size(otto,1));

                for x=1:size(otto)
                    set(otto(x),'MarkerEdgeColor',ottocal(x,:));
                    set(otto(x),'MarkerFaceColor',ottocal(x,:));
                end

                plot(trial_normalized,'k','LineWidth',2);
                plot(hem2trial_normalized, 'r','LineWidth',2);

                hold off;
                set(otto,'buttonDownFcn',@TC_Click_Trialselect_Callback);
                axis tight;

            else

                temp = find(tcm.session_ons_mean{1} == max(tcm.session_ons_mean{1}));
                otto = errorbar(tcm.trial_normalized{1}-100,tcm.trial_se_normalized{1});
                axis tight;
                line([temp,temp],[-10,10],'LineWidth',25,'Color',[1 0.5 0.5])

                otto = errorbar(tcm.trial_normalized{1}-100,tcm.trial_se_normalized{1},'k','LineWidth',2);
                xlim([0 length(tcm.trial_normalized{1})+1]);

                xlabel('Volume in trial');

            end

            psc_string = sprintf('%1.2f',Tempvar.psc);
            snr_string = sprintf('%1.2f',snr_trial);
            corr_string = sprintf('%1.2f',corr_trial);

            text(0.02,0.8,['Trial Corr = ' corr_string ],'Interpreter','none','FontSize',8,'Units','normalized','Color','Black','BackgroundColor','White');
            text(0.02,0.08,['signal change = ' psc_string ' %'],'Interpreter','none','FontSize',8,'Units','normalized','Color','Black','BackgroundColor','White');
            sandbox.textswitch = text(0.02,0.9,['Trial SNR = ' snr_string ],'Interpreter','none','FontSize',8,'Units','normalized','Color','Black','BackgroundColor','White');

            set(sandbox.textswitch,'buttonDownFcn',@TC_Display_Toggle_Callback);

            if analtrial ==1;
                sandbox.hemodynamics = hemodynamics_spare;
                sandbox.intrial = intrial_spare;
                idat = idat_spare;
            end

            Tempvar.feedb = ['Cluster size = ' num2str(cluster_size)];
            disp (['Cluster size = ' num2str(cluster_size)]);
            clear('tcm');
            refresh_feedb;
        end
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Click on the timecourse and select/unselect the corresponding trial
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function TC_Click_Trialselect_Callback(source,eventdata)

        franz = gco; % Get the lineseries object handle

        if strcmp(get(franz,'Selected'),'off')==1;

            set(franz,'Selected','on');
            seltr = find(franz==otto);
            disp (['Selected trial: ' num2str(seltr)]);
            temp_trial = find(sandbox.intrial(:,3)==seltr);
            sandbox.intrial(temp_trial,1)= 3;
            set(gcf,'CurrentAxes',hd); hold on;
            intrial2 = sandbox.intrial;
            intrial2(intrial2==0)=NaN;
            trialsel = bar (intrial2(:,1)*3+5,'y','EdgeColor','y');
            set(trialsel,'buttonDownFcn',@Click_Trialselect_Callback);
            hold off;
            refresh_image;

        else

            set(franz,'Selected','off');
            seltr = find(franz==otto);
            disp (['Deselected trial: ' num2str(seltr)]);
            temp_trial = find(sandbox.intrial(:,3)==seltr);
            sandbox.intrial(temp_trial,1)= 1;
            refresh_image;

        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Click on the SNR text to toggle between display modes
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function TC_Display_Toggle_Callback(source,eventdata)
        if ~isfield(sandbox,'plotall') || sandbox.plotall ==0;
            sandbox.plotall = 1;
        else
            sandbox.plotall =0;
        end
        Tempvar.imagechange =1;
        Trial_Overlay_Callback;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Call the brain extraction toolbox for further processing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Extract_Toolbox_Callback(source,eventdata)
        corint_idat = brainmask_toolbox(idat_spare,sandbox);
        refresh_map;
        idat = corint_idat;
        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end

        cla;
        refresh_volinfo;
        refresh_map;
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Normalize intensity over single volume (only x and y direction)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Correct_Intensity_Gradient_Callback(source,eventdata)

        %corint_idat = fsb_correct_intensity3(idat_spare,3); % New
        %experimental normalization 2
        %corint_idat = fsb_correct_intensity2(idat_spare,10,10); % New
        %experimental normalization 1
        corint_idat = fsb_correct_intensity(idat_spare,10,10); % Old version
        refresh_map;
        idat = corint_idat;
        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end

        cla;
        refresh_volinfo;
        refresh_map;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Apply a lens filter to the brain
% This is useful if brain extraction is not working
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Brain_Lens_Filter_Callback(source,eventdata)
        tic;
        Tempvar.lensfilt_param =get(h_lensfilt_slid,'Value');
        Tempvar.feedb = 'Filtering ongoing...';
        refresh_feedb;
        imean = fsb_calculate_mean(idat);
        [sandbox.wimage imean] = fsb_lens_brain (imean,Tempvar.lensfilt_param,Tempvar.slice_n);
        Tempvar.feedb =['Lensfilter threshold ' num2str(Tempvar.lensfilt_param) ];
        idat = imean;
        Tempvar.slice_n(4) = 1;
        refresh_figure;
        refresh_data;
        refresh_image;
        idat = idat_spare;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Realign brains
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Realign_Callback(source,eventdata)
        tic;

        str = get(source, 'String');
        val = get(source,'Value');

        if ~iscell(val)

            switch str;

                case 'MotCor' % Display CoM shift data
                    Tempvar.ctr_disp2 =0;

                case 'Toolbox' % Start up realignment toolbox
                    idat = fsb_realign(idat,sandbox);
            end

        else
            % Call function to run uncommon methods
            idat = fsb_realign_select(str,val,idat);

        end

        idat_spare = idat;
        dim = size(idat);
        t = toc;
        Tempvar.feedb =['Realignment done in ' num2str(t) ' seconds'];

        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end

        set(0,'CurrentFigure',sandbox_GUI);
        refresh_map;
        refresh_volinfo;
        refresh_all;
        Volume_Shift_COM_Callback;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Spatial smoothing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Filter_Idat_Spatial_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        tic;
        s_idat = fsb_ssmooth_idat(idat,Tempvar.ssmooth_param,1);
        idat = s_idat;
        cla;
        t = toc;
        Tempvar.feedb =['Spatial smoothing done in ' num2str(t) ' seconds'];
        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end
        idat_spare = idat;
        refresh_map;
        refresh_volinfo;
        refresh_all;
        Tempvar.imagechange = 0;
    end

    function Filter_Idat_Spatial_Slider_Callback(source,eventdata)
        Tempvar.ssmooth_vox = round((get(h_ssmooth_slider,'Value')));
        Tempvar.ssmooth_param = Tempvar.ssmooth_vox/2.355; % Formula for calculating the smoothing kernel as mentioned in Wikipedia(fwhm)
        Tempvar.feedb =['Spatial smoothing by ' num2str(Tempvar.ssmooth_vox) ' voxels'];
        disp(['Smoothing kernel Gaussian with ' num2str(Tempvar.ssmooth_vox) ' voxels extent'])
        refresh_image;
        refresh_feedb;
        refresh_text;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Temporal smoothing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Smooth_Idat_Temporal_Button_Callback(source,eventdata)
        tic;
        idat_spare = idat;
        h = waitbar(0,'Smoothing temporally...');
        waitbar(2/10)
        idat =permute(idat,[4 1 2 3]);
        waitbar(3/10)
        idat = single(idat(:));
        idat = smooth(idat,Tempvar.t_param);
        waitbar(5/10)
        idat =reshape(idat,dim(4),dim(1),dim(2),dim(3));
        waitbar(8/10)
        idat = permute(idat,[2 3 4 1]);
        waitbar(9/10)

        close(h) ;

        cla;

        t = toc;
        Tempvar.feedb =['Temporal smoothing done in ' num2str(t) ' seconds'];

        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end

        refresh_map;
        refresh_volinfo;
        refresh_all;
    end

    function Smooth_Idat_Temporal_Slider_Callback(source,eventdata)
        Tempvar.t_param = floor(get(h_tsmooth_slider,'Value'));
        Tempvar.feedb =['Temporal smoothing parameter ' num2str(Tempvar.t_param)];
        refresh_feedb;
        refresh_text;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Set low Filter Frequency
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Filter_low_Frequency_Callback(source,eventdata)

        filt_freq = get(h_filter_low,'String');
        set(h_filter_low,'String',num2str(filt_freq));
        sandbox.filt_freq_low = 1/str2double(filt_freq);

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Set high Filter Frequency
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Filter_high_Frequency_Callback(source,eventdata)

        filt_freq = get(h_filter_high,'String');
        set(h_filter_high,'String',num2str(filt_freq));
        sandbox.filt_freq_high = 1/str2double(filt_freq);

    end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Temporally filtering the data
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Filter_Temporal_Callback(source,eventdata)
        str = get(source, 'String');
        val = get(source,'Value');

        switch str{val};

            case 'Spectrum'
                sandbox.filttype = 'Spectrum';
                sandbox.Tempvar.roi = Tempvar.roi;
                sandbox.Tempvar.slice_n(1) = Tempvar.slice_n(1);
                sandbox.Tempvar.slice_n(2) = Tempvar.slice_n(2);
                sandbox.Tempvar.slice_n(3) = Tempvar.slice_n(3);
                sandbox.filt_freq = str2double(get(h_filter_high,'String'));
            case 'Low-Pass Filter'
                sandbox.filttype = 'Low';
                sandbox.filt_freq = str2double(get(h_filter_low,'String'));

            case 'High-Pass Filter'
                sandbox.filttype = 'High';
                sandbox.filt_freq = str2double(get(h_filter_high,'String'));

            case 'Band-Pass Filter'
                sandbox.filttype = 'Stop';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;

            case 'Band-Pass Filter run by run'
                sandbox.filttype = 'Stop';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;


            case 'Remove discontinuities';
                idat = fsb_remove_disc(idat,sandbox);
                refresh_all;
                return
            case 'Interpolate out of trial periods and filter data'
                disp('Interpolating artefacts run by run and applying filter...');
                sandbox.filttype = 'Stop';
                sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                    str2double(get(h_filter_low,'String'))] ;
                idat = fsb_interpolate_filter_scans(idat,sandbox);
                %                filtopts.discard = 0;
                %                filtopts.invert = 0;
                %                filtopts.filter = 0;
                %                filtopts.mintrial = 5;
                %                idat = fsb_interpolate_filter(idat,sandbox,filtopts);
                refresh_all
                return
                
                %             case 'SPM filter' % Currently not implemented
                %                 sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                %                     str2double(get(h_filter_low,'String'))] ;
                %                 idat = fsb_spm_filter(idat,sandbox);
                %                 refresh_all;
                %                 return
                %
                %             case 'Interpolate Artefacts and filter data'
                %                 disp('Interpolating artefacts run by run and applying filter...');
                %                 sandbox.filttype = 'Stop';
                %                 sandbox.filt_freq = [str2double(get(h_filter_high,'String'))...
                %                     str2double(get(h_filter_low,'String'))] ;
                %                 [idat sandbox] = fsb_interpolate_filter_artefacts(idat,sandbox);
                %                 refresh_all
                %                 return

        end

        idat = fsb_filter_temporal_scan(idat,sandbox,Tempvar);
        refresh_map;
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Display slices in another window
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function  Slice_display_Radiobutton_Callback(source,eventdata)

        if (get(h_showslices,'Value') == get(h_showslices,'Max'))
            param.dispslices = 1;

            if isempty(Map.maptype)
                Map.maptype = 'Plot slices';
            end
        else
            param.dispslices = 0;
            close(103);
            %Map.maptype = [];
        end
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate maps for whole brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    function  Map_calculate_volume_Callback(source,eventdata)
        tic;
        if ~isempty(sandbox.hemodynamics) && size(sandbox.hemodynamics,1) ==dim(4)
            Map.corparam =get(h_calcorrelation2_slid,'Value');
            idat2=idat;
            all_mean = mean(imean(:));
            idat2(idat2<all_mean/2) = 0;
            sandbox = fsb_calculate_full_map(sandbox,idat2,param.hrf_pred,Map.statsel);
            if Tempvar.clustparam>1;
                if Tempvar.sparecor ==0;

                    r_cor_red(:,:,:,param.hrf_pred)= fsb_threshold_cluster(r_cor_red(:,:,:,param.hrf_pred),Tempvar.clustparam);
                    r_cor_neg(:,:,:,param.hrf_pred) = fsb_threshold_cluster(r_cor_neg(:,:,:,param.hrf_pred),Tempvar.clustparam);

                else

                    r_cor_red_spare = r_cor_red;
                    r_cor_neg_spare = r_cor_neg;
                    Tempvar.sparecor =1;

                end
            end

            r_cor_red = sandbox.stats.r_cor>Map.corparam;
            r_cor_neg =sandbox.stats.r_cor <(Map.corparam*(-1));
            t = toc;
            Tempvar.feedb =[' done in ' num2str(t) ' seconds'];
        else
            Tempvar.feedb = 'No hemodynamics data loaded';
        end
        Tempvar.mapcalc = 1;
        Tempvar.imagechange=0;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Get statistics setup
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Map_display_Button_Callback(source,eventdata)

        if (get(h_mapdisp,'Value') == get(h_mapdisp,'Max'))
            Tempvar.mapdisp = 1;
            Map_threshold_slider_Callback;
        else
            Tempvar.mapdisp = 0;
        end
        Tempvar.imagechange = 1;
        refresh_all;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate correlation map for displayed slices only
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


    function  Map_threshold_edit_Callback(source,eventdata)

        if size(sandbox.hemodynamics,1) ==dim(4) && Tempvar.mapdisp ==1;

            Map.corparam = str2num(get(h_probability_button,'string'));
            Map.probparam = fsb_calculate_prob(dim(4),Map.corparam);
            probparams = sprintf('%1.5f',Map.probparam);
            corparams = sprintf('%1.4f',Map.corparam);
            probstring = ['Probability = ' probparams];
            set(h_calcorrelation2_slid,'Value');
            set(h_probability_slid,'String',probstring);
            disp([' Correlation threshold ' num2str(Map.corparam) ]);
        else
            Tempvar.feedb = 'No hemodynamics data loaded';
            disp('No hemodynamics data loaded');
        end

        refresh_all;

    end

    function  Map_threshold_slider_Callback(source,eventdata)

        if size(sandbox.hemodynamics,1) ==dim(4) && Tempvar.mapdisp ==1;
            Map.corparam =get(h_calcorrelation2_slid,'Value');
            Map.probparam = fsb_calculate_prob(dim(4),Map.corparam);
            probparams = sprintf('%1.5f',Map.probparam);
            corparams = sprintf('%1.4f',Map.corparam);
            probstring = ['Probability = ' probparams];
            set(h_probability_button,'String',corparams);
            set(h_probability_slid,'String',probstring);
            disp([' Correlation threshold ' num2str(Map.corparam) ]);
        else
            Tempvar.feedb = 'No hemodynamics data loaded';
            disp('No hemodynamics data loaded');
        end

        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate Sandbox statistical map for displayed slices only
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Map_calculate_red
        tic;
        idat2=idat;
        all_mean = mean(imean(:));
        idat2(idat2<all_mean/2) = 0;
        sandbox = fsb_calculate_red_map(sandbox,idat2,Map.statsel,Tempvar.slice_l,param.hrf_pred);
        t=  toc; disp(['Map calculation ' num2str(t) ' seconds']);
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Analyze only selected trials
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Trial_analyze_selected_Callback(source,eventdata)

        if  isempty(find (sandbox.intrial(:,1)==3));
            disp('No trial selected, aborting')
            return
        end

        if (get(h_trialselrb,'Value') == get(h_trialselrb,'Max'))

            analtrial = 1;
            temp = find (sandbox.intrial(:,1)==3);
            idat_red = idat(:,:,:,temp);

            r_cor_spare =sandbox.stats.r_cor;
            r_cor_red_spare = r_cor_red;
            r_cor_neg_spare =r_cor_neg;
            hemodynamics_spare = sandbox.hemodynamics;
            hemodynamics_red = sandbox.hemodynamics(temp,:);

            intrial_spare = sandbox.intrial;
            onset_spare = sandbox.onset;

            sandbox.stats.r_cor = zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
            r_cor_red = zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
            r_cor_neg =zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
        else
            analtrial = 0;
            sandbox.stats.r_cor =r_cor_spare;
            r_cor_red = r_cor_red_spare;
            r_cor_neg =r_cor_neg_spare;
            sandbox.hemodynamics=hemodynamics_spare;
        end

        refresh_map;
        refresh_all;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate correlation map for selected trials only
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function  Map_calculate_red_trial

        selection = sandbox.intrial(:,1)>1;
        sandbox = fsb_calculate_red_map(sandbox,idat_red,Map.statsel,Tempvar.slice_l,param.hrf_pred,selection);

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Get clustering setup
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Map_cluster_size_Slider_Callback(source,eventdata)
        if size(sandbox.hemodynamics,1) ==dim(4) && Tempvar.mapdisp ==1;
            Tempvar.clustparam =get(h_clustering_slid,'Value');
            Tempvar.feedb =[' Clustering threshold ' num2str(Tempvar.clustparam) ];

        else

            Tempvar.feedb = 'No hemodynamics data loaded';

        end
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% shift hemodynamics back and forth with edit field
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HRF_edit_Callback(source,eventdata)

        if size(sandbox.hemodynamics,1) ~=dim(4);
            Keep_Modifications
        end
        try
            inparg = struct;
            Tempvar.imagechange = 1;
            intrial_spare(:,20:end)=[];
            onset_spare(:,20:end)=[];
            hemodynamics_spare(:,20:end)=[];
            if param.hrf_pred>20
                param.hrf_pred =1;
            end

            sandbox.intrial = intrial_spare;
            sandbox.onset = onset_spare;
            sandbox.hemodynamics = hemodynamics_spare;

            Tempvar.shiftparam =str2num(get(h_shifthrfedit,'string'));
            Tempvar.feedb =['Shiftparameter ' num2str(Tempvar.shiftparam) ' Volume' ];
            disp(['Shiftparameter ' num2str(Tempvar.shiftparam) ' Volume' ]);

            try
                for x=1:size(sandbox.volmat.evt_onsets,2);
                    inparg.evt_onsets{x,1} = sandbox.volmat.evt_onsets{x}+Tempvar.shiftparam;
                end

                inparg.evt_names = sandbox.volmat.evt_names;
                inparg.evt_dur = sandbox.volmat.evt_dur;
                inparg.TR = sandbox.TR;
                inparg.volumes = dim(4);
                inparg.deriv = 0;
                [U,bf,p,X] = fsb_conv_onsets(inparg);
                sandbox.hemodynamics = X;
            catch % if there is no onset information present
                sandbox.hemodynamics = circshift(sandbox.hemodynamics,[round(Tempvar.shiftparam) 0]);
            end
            disp('Onset convolution done');
            sandbox.intrial = circshift(sandbox.intrial,[round(Tempvar.shiftparam) 0]);
            disp('Intrial shifting done');
            sandbox.onset = circshift(sandbox.onset,[round(Tempvar.shiftparam) 0]);
            disp('Onset shifting done');
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            Tempvar.feedb = 'Problems detected, operation aborting';
            disp('Problems detected, operation aborting');
        end

        refresh_map;
        refresh_all;
        Tempvar.imagechange=0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% shift hemodynamics back and forth
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HRF_shift_Slider_Callback(source,eventdata)

        if size(sandbox.hemodynamics,1) ~=dim(4);
            Keep_Modifications
        end
        try
            inparg = struct;
            Tempvar.imagechange = 1;
            intrial_spare(:,20:end)=[];
            onset_spare(:,20:end)=[];
            hemodynamics_spare(:,20:end)=[];
            if param.hrf_pred>20
                param.hrf_pred =1;
            end
            sandbox.intrial = intrial_spare;
            sandbox.onset = onset_spare;
            sandbox.hemodynamics = hemodynamics_spare;


            Tempvar.shiftparam =(get(h_shifthrfslider,'Value'));
            Tempvar.feedb =['Shiftparameter ' num2str(Tempvar.shiftparam) ' Volume' ];
            disp(['Shiftparameter ' num2str(Tempvar.shiftparam) ' Volume' ]);

            for x=1:size(sandbox.volmat.evt_onsets,2);
                inparg.evt_onsets{x,1} = sandbox.volmat.evt_onsets{x}+Tempvar.shiftparam;
            end
            inparg.evt_names = sandbox.volmat.evt_names;
            inparg.evt_dur = sandbox.volmat.evt_dur;
            inparg.TR = sandbox.TR;
            inparg.volumes = dim(4);
            inparg.deriv = 0;
            [U,bf,p,X] = fsb_conv_onsets(inparg);
            sandbox.hemodynamics = X;
            disp('Onset convolution done');
            sandbox.intrial = circshift(sandbox.intrial,[round(Tempvar.shiftparam) 0]);
            disp('Intrial shifting done');
            sandbox.onset = circshift(sandbox.onset,[round(Tempvar.shiftparam) 0]);
            disp('Onset shifting done');
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            Tempvar.feedb = 'Problems detected, operation aborting';
            disp('Problems detected, operation aborting');
        end

        refresh_map;
        refresh_all;
        Tempvar.imagechange=0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% reset hemodynamic shift
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HRF_shift_Button_Callback(source,eventdata)
        Tempvar.imagechange =1;
        intrial_spare(:,20:end)=[];
        onset_spare(:,20:end)=[];
        hemodynamics_spare(:,20:end)=[];
        if param.hrf_pred>20
            param.hrf_pred =1;
        end
        sandbox.intrial = intrial_spare;
        sandbox.onset = onset_spare;
        sandbox.hemodynamics = hemodynamics_spare;

        Tempvar.shiftparam =0;
        set(h_shifthrfslider,'Value',0);
        try
            sandbox.volmat.evt_onsets = sandbox.origdata.volmat.evt_onsets;
            sandbox.intrial = sandbox.origdata.intrial;
            sandbox.hemodynamics = sandbox.origdata.hemodynamics;
            sandbox.timemat.evt_onsets = sandbox.origdata.timemat.evt_onsets;
            sandbox.onset = sandbox.origdata.onset;
            inparg.evt_onsets = sandbox.volmat.evt_onsets;

            inparg.evt_dur = sandbox.volmat.evt_dur;
            %inparg.evt_names = sandbox.volmat.evt_names;
            inparg.TR = sandbox.TR;
            inparg.volumes = dim(4);
            inparg.deriv = 0;
            Tempvar.feedb =['Shiftparameter reset to ' num2str(Tempvar.shiftparam) ' Volumes' ];
            disp(['Shiftparameter reset to ' num2str(Tempvar.shiftparam) ' Volumes' ]);
            [U,bf,p,X] = fsb_conv_onsets(inparg);
            sandbox.hemodynamics = X;

            refresh_map;
            refresh_all;
            Tempvar.imagechange=0;
        catch
            disp('No prior optimization found')
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% look for the optimal shift of your predictor
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HRF_auto_adjust_Callback(source, eventdata)

        Tempvar.imagechange = 1;
        if size(sandbox.hemodynamics,1) ~=dim(4);
            Keep_Modifications
        end
        param.Tempvar.roi = Tempvar.roi;
        param.Tempvar.slice_n = Tempvar.slice_n;
        try
            sandbox = fsb_HRF_auto_adjust(idat,sandbox,param);
        catch
            errmsg = lasterror;
            fsb_error_handler(errmsg);
            disp('Hemodynamic predictor shifting failed, please check predictors')
        end

        refresh_map;
        refresh_all;
        Tempvar.imagechange=0;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% smooth hemodynamics predictors
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function HRF_smooth_Slider_Callback(source,eventdata)

        if size(sandbox.hemodynamics,1) ==dim(4) && Tempvar.mapdisp ==1;
            Tempvar.smoothparam =round(get(h_smoothhrfslider,'Value'));
            Tempvar.feedb =['HRF smoothing parameter ' num2str(Tempvar.smoothparam) ' Volume         ' ];

        else

            Tempvar.feedb = 'No hemodynamics data loaded';

        end

        refresh_feedb;

    end

    function HRF_smooth_Button_Callback(source,eventdata)
        Tempvar.imagechange = 1;
        Tempvar.smoothed =1;
        if size(sandbox.hemodynamics,1) ==dim(4) && Tempvar.mapdisp ==1;
            Tempvar.smoothed = Tempvar.smoothed +Tempvar.smoothparam;

            if sum(hemodynamics_spare) ==0;
                hemodynamics_spare = sandbox.hemodynamics;
            end

            sandbox.hemodynamics = smooth(sandbox.hemodynamics,Tempvar.smoothparam);
            sandbox.hemodynamics = reshape(sandbox.hemodynamics,dim(4),[]);
            Tempvar.feedb =['Tempvar.smoothed by ' num2str(Tempvar.smoothed) ' Volumes' ];
            if analtrial ==1;
                idat_red = idat(:,:,:,temp);
                hemodynamics_red = hemodynamics_spare(temp,:);

            end

        else
            Tempvar.feedb = 'No hemodynamics data loaded';

        end
        refresh_map;
        refresh_all;
        Tempvar.imagechange = 0;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Undo
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Undo_Button_Callback(source,eventdata)
        sandbox.hemodynamics = hemodynamics_spare;
        idat = idat_spare;
        sandbox.onset = onset_spare;
        sandbox.intrial = intrial_spare;

        if analtrial ==1;
            idat_red = idat(:,:,:,temp);
            hemodynamics_red = hemodynamics_spare(temp,:);
        end

        refresh_volinfo;
        refresh_map;
        refresh_all;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Read the helpfile
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Help_Button_Callback(source, eventdata)
        web sandbox_help.html -helpbrowser
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% center of mass calculation -  copyright
% George A. Keliris 2006
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function Volume_Shift_COM_Callback(source,eventdata)
        tic;
        if Tempvar.ctr_disp2 ==0;

            [dctr_mass_r, dctr_mass_xyz]=fsb_ctrm_3D_sandbox(idat,1,0,sandbox);
            ctr_disp = dctr_mass_r*(dim(1)/3);
            ctr_mass =[];
            Tempvar.ctr_disp2 =1;
            ctr_mass(:,1)=dctr_mass_xyz(:,2);
            figure(104); set(gcf,'Name','Movement parameters');
            clf;
            lw = 150/ length(sandbox.intrial);
            lw =lw*3;
            %lw =lw*5;
            temp = find(sandbox.intrial(:,1) >= 1);
            o4 = line([temp,temp],[-6,14],'LineWidth',lw,'Color',[1 1 0.5]);
            set(o4,'buttonDownFcn',@Click_Time_Callback);
            hold on;

            o1 = plot(dctr_mass_xyz(:,2),'LineWidth',2);
            set(o1,'buttonDownFcn',@Click_Time_Callback);

            if isfield(sandbox,'sensdata') && ~isempty(sandbox.sensdata);

                o2 = plot (sandbox.sensdata(:,3),'g');
                set(o2,'buttonDownFcn',@Click_Time_Callback);
                o3 = plot (sandbox.sensdata(:,4),'r');

                set(o3,'buttonDownFcn',@Click_Time_Callback);
                axis tight;
                legend([o1 o2 o3 o4(1)],'COM','Sensor A', 'Sensor B','Trial periods', 'location', 'NorthWest');
            else
                axis tight;
                legend([o1 o4(1)],'COM','Trial periods', 'location', 'NorthWest');

            end

            xlabel('time -> [img]');
            ylabel('center-of-mass shift -> [mm]');
            Title('Estimated center-of-mass shift (A-P)');
            disp (['Max. negative translation: ' num2str(min(ctr_mass)) ' mm']);
            disp (['Max. positive translation: ' num2str(max(ctr_mass)) ' mm']);
        else
            Tempvar.ctr_disp2=0;
        end
        refresh_image;
        t=toc; disp(['Center of Mass calculation done in ' num2str(t) ' s']);
    end


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% refresh the data
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_data

        if Tempvar.hoverlay ==1;
            slicex_data = squeeze(hidat(Tempvar.slice_n(1),:,:));
            slicex_data = slicex_data';
            slicey_data = squeeze(hidat(:,Tempvar.slice_n(2),:));
            slicey_data = slicey_data';
            slicez_data = squeeze(hidat(:,:,Tempvar.slice_n(3)));
            Tempvar.slice_h = Tempvar.slice_n;
            Tempvar.slice_l(1:3) = round(Tempvar.slice_n(1:3)./Tempvar.dimratio);
            slicet_data = squeeze(idat(Tempvar.slice_l(1),:,Tempvar.slice_l(3),:));
        else
            slicex_data = squeeze(idat(Tempvar.slice_n(1),:,:,Tempvar.slice_n(4)));
            slicex_data = slicex_data';
            slicey_data = squeeze(idat(:,Tempvar.slice_n(2),:,Tempvar.slice_n(4)));
            slicey_data = slicey_data';
            slicez_data = squeeze(idat(:,:,Tempvar.slice_n(3),Tempvar.slice_n(4)));
            slicet_data = squeeze(idat(Tempvar.slice_n(1),:,Tempvar.slice_n(3),:));
            Tempvar.slice_l = Tempvar.slice_n;
            Tempvar.slice_h(1:3) = round(Tempvar.slice_n(1:3).*Tempvar.dimratio);

        end

    end

    function refresh_buttons

        set(h_showslices,'Value',0);
        set(h_overlaytrial,'Value',0);
        set(h_overlayclus,'Value',0);
        try
            set(h_highres,'Value',0);
        catch
        end
        set(h_mapdisp,'Value',0);
        set(h_trialselrb,'Value',0);
        set(h_showslices,'Value',0);

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% refresh feedback
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_feedb
        if Tempvar.feedb >1;
            set(gcf,'CurrentAxes',feedb_win);
            cla;
            try
                text(0.01,0.01,Tempvar.feedb,'Interpreter','none','FontSize',12,'Units','normalized','Color','Black','BackgroundColor','White');
            catch
                disp(['Tempvar.feedb='  num2str( Tempvar.feedb)]);
            end
            drawnow expose;
            Tempvar.feedb = 0;
        end
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% refresh text
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_text

        set(0,'CurrentFigure',sandbox_GUI);
        set (gcf,'CurrentAxes',info_win); axis off; axis ij;cla;
        text(0,0.7,['x= ' num2str(sandbox.hdr.dime.pixdim(2)) 'mm'],'FontSize',8,'Units','normalized','Color','Black');
        text(0,0.6,['y= ' num2str(sandbox.hdr.dime.pixdim(3)) 'mm'],'FontSize',8,'Units','normalized','Color','Black');
        text(0,0.5,['z= ' num2str(sandbox.hdr.dime.pixdim(4)) 'mm'],'FontSize',8,'Units','normalized','Color','Black');

        try
            set (gcf,'CurrentAxes',beta_win);
            cla;
            bar('v6',best,'r');axis off; axis xy; hold on;
            my_errorbar(1:length(best)+.3,best,best-bint(:,1),bint(:,2)-best);
            hold off;
        catch

        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Refresh the image
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_image

        set(0,'CurrentFigure',sandbox_GUI);
        %drawnow;
        maskit =1;
        set(gcf,'CurrentAxes',ha)
        image_x=imagesc(slicex_data);
        axis xy; axis off;
        set(image_x,'buttonDownFcn',@Click_x_Callback);
        text(0.01,0.1,num2str(Tempvar.slice_n(1)),'FontSize',10,'Units',...
            'normalized','Color','White','FontWeight','Bold');
        Tempvar = fsb_ROI_mask(idat,Tempvar);

        if  Tempvar.mapdisp ==1

            if size(sandbox.stats.r_cor,4)<param.hrf_pred;
                sandbox.stats.r_cor(:,:,:,param.hrf_pred) = 0;
            end

            if analtrial ==0
                if Tempvar.imagechange == 1 && Tempvar.mapcalc == 0;
                    Map_calculate_red;
                end
            else
                Map_calculate_red_trial;
            end

            switch Map.statsel;

                case 'corrc'
                    r_cor_red = sandbox.stats.r_cor >Map.corparam;
                    r_cor_neg = sandbox.stats.r_cor <(Map.corparam*(-1));
                    r_mask = zeros(size(r_cor_red));
                    r_mask(r_cor_red) =1;
                    r_mask(r_cor_neg) =1;
                    dummy = sandbox.stats.r_cor.*r_mask;
                    Map.rc =1;

                case 'glm';

                    r_cor_red = sandbox.stats.glm_prob <Map.probparam;
                    r_cor_neg =sandbox.stats.glm_prob  >1;
                    dummy = r_cor_red.*sandbox.stats.glm_prob;
                    Map.rc=0;

                case 'robreg'
                    r_cor_red = sandbox.stats.rob_prob <Map.probparam;
                    r_cor_neg =sandbox.stats.rob_prob  >1;
                    dummy = r_cor_red.*sandbox.stats.rob_prob;
                    Map.rc=0;
                case 'ICA'
                    r_cor_red = dummy >0;
                    r_cor_neg =dummy<1;
                    Map.rc=99;

            end
            
            hold on;

            if Tempvar.clustparam>1;

                r_cor_red(:,:,:,param.hrf_pred)= fsb_threshold_cluster(r_cor_red(:,:,:,param.hrf_pred),Tempvar.clustparam);
                r_cor_neg(:,:,:,param.hrf_pred) = fsb_threshold_cluster(r_cor_neg(:,:,:,param.hrf_pred),Tempvar.clustparam);

                if param.hrf_pred==0;
                    param.hrf_pred =1;
                end

                r_mask = zeros(size(r_cor_red(:,:,:,size(param.hrf_pred,1))));
                r_mask(r_cor_red(:,:,:,param.hrf_pred)>0) =1;
                r_mask(r_cor_neg(:,:,:,param.hrf_pred)>0) =1;
                dummy(:,:,:,param.hrf_pred) = dummy(:,:,:,param.hrf_pred).*r_mask;

            end

            Map.z = (squeeze(dummy(Tempvar.slice_l(1),:,:,param.hrf_pred))');

            Map = fsb_colmap(Map);
            
            Map.slice_data = slicex_data;

            if Tempvar.hoverlay == 1;
                Map = fsb_map_overlay(Map);
            end

            plot_x1 = image(Map.zm);
            set (plot_x1,'AlphaData',Map.tzi);
            set(plot_x1,'buttonDownFcn',@Click_x_Callback);

            hold off;

        end

        line([Tempvar.slice_n(2),Tempvar.slice_n(2)],[0,512],'LineWidth',1,'Color','m')
        line([512,0],[Tempvar.slice_n(3),Tempvar.slice_n(3)],'LineWidth',1,'Color','m')
        rectangle('Position',[Tempvar.slice_n(2)-Tempvar.roi,Tempvar.slice_n(3)-Tempvar.roi,Tempvar.roi*2,Tempvar.roi*2],'EdgeColor','g');
        %rectangle('Position',[Tempvar.roiind(1,2),Tempvar.roiind(1,3),Tempvar.roi*2,Tempvar.roi*2],'EdgeColor','g');

        set(gcf,'CurrentAxes',hb)
        image_y=imagesc(slicey_data);
        axis xy; axis off;
        text(0.01,0.1,num2str(Tempvar.slice_n(2)),'FontSize',10,'Units','normalized','Color','White','FontWeight','Bold');
        set(image_y,'buttonDownFcn',@Click_y_Callback);

        if Tempvar.mapdisp ==1;
            hold on;

            Map.z = (squeeze(dummy(:,Tempvar.slice_l(2),:,param.hrf_pred)))';

            Map = fsb_colmap(Map);
            Map.slice_data = slicey_data;

            if Tempvar.hoverlay == 1;
                Map = fsb_map_overlay(Map);
            end

            plot_y1 = image(Map.zm);
            set (plot_y1,'AlphaData',Map.tzi);
            set(plot_y1,'buttonDownFcn',@Click_y_Callback);
            hold off;
        end

        line([Tempvar.slice_n(1),Tempvar.slice_n(1)],[0,512],'LineWidth',1,'Color','m')
        line([512,0],[Tempvar.slice_n(3),Tempvar.slice_n(3)],'LineWidth',1,'Color','m')
        rectangle('Position',[Tempvar.slice_n(1)-Tempvar.roi,Tempvar.slice_n(3)-Tempvar.roi,Tempvar.roi*2,Tempvar.roi*2],'EdgeColor','g');

        set(gcf,'CurrentAxes',hc)
        image_z=imagesc(slicez_data);
        axis xy; axis off;
        text(0.01,0.05,num2str(Tempvar.slice_n(3)),'FontSize',10,'Units','normalized','Color','White','FontWeight','Bold');


        set(image_z,'buttonDownFcn',@Click_z_Callback);

        if Tempvar.mapdisp ==1;

            hold on;
            Map.z = squeeze(dummy(:,:,Tempvar.slice_n(3),param.hrf_pred));
            Map = fsb_colmap(Map);
            Map.slice_data = slicez_data;

            if Tempvar.hoverlay == 1;
                Map = fsb_map_overlay(Map);
            end

            plot_z1 = image(Map.zm);
            set (plot_z1,'AlphaData',Map.tzi);
            set(plot_z1,'buttonDownFcn',@Click_z_Callback);
            hold off;

        end
        voxel_value = idat(Tempvar.slice_l(1),Tempvar.slice_l(2),Tempvar.slice_l(3),Tempvar.slice_n(4));
        voxel_value_string = sprintf('%1.2f',voxel_value);

        if Tempvar.hoverlay ==1;
            text(0.29,1.00,['Voxel Value: ' num2str(voxel_value_string)],...
                'FontSize',9,'Units','normalized','Color','Magenta','BackgroundColor','Black');
        else
            text(0.1,1.0,['Voxel Value: ' num2str(voxel_value_string)],...
                'FontSize',9,'Units','normalized','Color','Magenta','BackgroundColor','Black');

            %roi_vox = mean(imean(Tempvar.roimask>0));
            idat_temp = idat(:,:,:,Tempvar.slice_n(4));
            roi_vox = mean(idat_temp(Tempvar.roimask>0));
            avg_roi = mean(roi_vox(:));
            avg_roi_string = sprintf('%1.2f',avg_roi);
            text(0.6,1.00,['ROI Mean: ' num2str(avg_roi_string)],'FontSize',9,'Units','normalized','Color','Green','BackgroundColor','Black');

        end

        line([Tempvar.slice_n(2),Tempvar.slice_n(2)],[1,512],'LineWidth',1,'Color','m')
        line([512,1],[Tempvar.slice_n(1),Tempvar.slice_n(1)],'LineWidth',1,'Color','m')
        rectangle('Position',[Tempvar.slice_n(2)-Tempvar.roi,Tempvar.slice_n(1)-Tempvar.roi,Tempvar.roi*2,Tempvar.roi*2],'EdgeColor','g');

        if Tempvar.mapdisp ==1;
            try
                set(gcf,'CurrentAxes',colmap); axis off;
                cm2 = flipud(Map.cm);
                cm2 = repmat(cm2,[1 1 10]);
                cm2 = permute(cm2,[1 3 2]);
                image(cm2);
                set(gca,'YTick',1:25:101);
                set(gca,'YTickLabel',[1 0.5  0 -0.5  -1],'FontSize',8);
                set(gca,'XTickLabel',' ');
            catch
                cla;
            end
        else
            set(gcf,'CurrentAxes',colmap); axis off; cla;
        end

        if size(dim,2)<4;

            hrf_avg_disp = Tempvar.roimean;
            avg_hrf_avg = mean(hrf_avg_disp);
            noise_std = std(single(imean(imean<mean(mean(mean(imean)))/4)));
            snr_roi = avg_hrf_avg/noise_std;
            snr_roi_string = sprintf('%1.2f',snr_roi);
            set(gcf,'CurrentAxes',hd);axis xy;drawnow expose;cla;
            text(0.0,1.05,['ROI SNR : ' snr_roi_string],'FontSize',9,'Units','normalized','Color','Black','BackgroundColor','White');

        else

            set(gcf,'CurrentAxes',hd);
            image_t=imagesc(slicet_data);
            axis xy;
            text(0.003,0.94,num2str(Tempvar.slice_n(4)),'FontSize',12,'Units','normalized','Color','White','FontWeight','Bold');
            set(image_t,'buttonDownFcn',@Click_Time_Callback);
            hold on;
            line([Tempvar.slice_n(4),Tempvar.slice_n(4)],[1,512],'LineWidth',1,'Color',[1 1 1])
            hrf_disp = single(squeeze(idat(Tempvar.slice_l(1),Tempvar.slice_l(2),Tempvar.slice_l(3),:)));
            min_hrf = min(hrf_disp);
            max_hrf = max(hrf_disp);
            avg_hrf = mean(hrf_disp);
            psc_hrf = (max_hrf-min_hrf)/avg_hrf*100;

            psc_avg_string = sprintf('%1.2f',Tempvar.psc_hrf_avg);
            text(0.6,1.05,['ROI change : ' psc_avg_string '%'],'FontSize',9,'Units','normalized','Color','Green','BackgroundColor','Black');

            % Noise area calculation based on image size
            noise_lim = round(size(idat)/10);
            noise_lim = noise_lim(1:2);
            noise_area = single(imean(1:noise_lim(1),1:noise_lim(2),:));

            % Noise area calculation based on image intensity
            %noise_area = single(imean(imean<mean(mean(mean(imean)))/4));

            noise_mean = mean(noise_area(:));
            %noise_std = std(noise_area(:));

            % SNR definition 1 according to ReviseMRI.com
            % SNR = avg(roi image)/ std(roi noise)
            % snr_roi = avg_hrf_avg/noise_std;

            % SNR definition 2 according to MRI - From Picture to Proton
            % McRobbie et al, Cambridge Univ. Press, 2003
            % SNR = avg(roi Image) / avg(roi Noise)
            snr_roi = Tempvar.avg_hrf_avg/noise_mean;

            snr_roi_string = sprintf('%1.2f',snr_roi);
            text(0.0,1.05,['ROI SNR : ' snr_roi_string],'FontSize',9,'Units','normalized','Color','Black');

            hrf_avg_disp2 = (Tempvar.max_hrf_avg-Tempvar.min_hrf_avg)/dim(2);

            %if hrf_avg_disp2 >0;
            hrf_avg_disp = (Tempvar.roimean-min(Tempvar.roimean))/hrf_avg_disp2;

            if min(hrf_avg_disp)<0
                hrf_avg_disp=hrf_avg_disp+abs(min(hrf_avg_disp));
                hrf_avg_disp2 = max(hrf_avg_disp)/dim(2);
                hrf_avg_disp = hrf_avg_disp/hrf_avg_disp2;
            end

            hrfdisp = plot(hrf_avg_disp,'LineWidth',2,'Color','g');
            set(hrfdisp,'buttonDownFcn',@Click_Time_Callback);
            %   end

            if min(hrf_disp)<0;
                hrf_disp=hrf_disp+abs(min(hrf_disp));
            end

            hrf_disp2 = max(hrf_disp)/dim(2);


            if hrf_disp2 >0;
                hrf_disp = (hrf_disp/hrf_disp2);
            end

            psc_string = sprintf('%1.2f',psc_hrf);

            text(0.4,1.05,['Voxel change : ' psc_string '%'],'FontSize',9,'Units','normalized','Color','Magenta','BackgroundColor','Black');
            text(0.8,1.05,['Volume change : ' psc_vol_string '%'],'FontSize',9,'Units','normalized','Color','Black');
            text(0.2,1.05,['SNR : ' num2str(Tempvar.snr) ],'FontSize',9,'Units','normalized','Color','Black');

            hrfdispmin=min(hrf_disp);
            hrfdispmax=max(hrf_disp);
            hrf_disp_single = (hrfdispmax-hrfdispmin)/dim(2);
            hrf_disp = (hrf_disp-hrfdispmin)/hrf_disp_single;


            hrfdisp2 = plot(hrf_disp,'LineWidth',1,'Color','m');
            set(hrfdisp2,'buttonDownFcn',@Click_Time_Callback);

            if Tempvar.ctr_disp2 >0;
                plot (ctr_disp, 'LineWidth', 2, 'Color',[0 0 1]);
            end

            if Tempvar.disp_disc == 1;
                plot (sandbox.discard,10,'-w+',...
                    'LineWidth',1,...
                    'MarkerEdgeColor','w',...
                    'MarkerFaceColor',[1 1 1],...
                    'MarkerSize',10)
            end

            if disp_onset == 1 && size(sandbox.hemodynamics,1)==dim(4) && max(hrf_disp)>0;
                hemmean = max(sandbox.hemodynamics(:,param.hrf_pred));
                bhem2 = (sandbox.hemodynamics(:,param.hrf_pred)/hemmean)*(dim(2)/4)+(dim(2)/4);

                hemdispmin=min(bhem2);
                hemdispmax=max(bhem2);
                hem_disp_single = (hemdispmax-hemdispmin)/dim(2);
                bhem2 = (bhem2-hemdispmin)/hem_disp_single;

                hemplot = plot(bhem2,...
                    '--ro','LineWidth',1,...
                    'MarkerEdgeColor','r',...
                    'MarkerFaceColor','k',...
                    'MarkerSize',3);hold off;
                set(hemplot,'buttonDownFcn',@Click_Time_Callback);

                if ~isempty(sandbox.intrial)
                    hold on;
                    intrial2 = sandbox.intrial;
                    intrial2(intrial2==0)=NaN;
                    trialsel = plot(intrial2(:,1),'y',...
                        'LineWidth',10);
                    set(trialsel,'buttonDownFcn',@Click_Time_Callback);
                    hold off;
                end

                if ~isempty(sandbox.onset);
                    hold on;
                    onset2 = sandbox.onset;
                    onset2(onset2==0)=NaN;
                    onsetsel = plot (onset2,'-b+',...
                        'LineWidth',1,...
                        'MarkerSize',10);
                    set(onsetsel,'buttonDownFcn',@Click_Time_Callback);
                    hold off;
                end

                if Tempvar.overlayt==1
                    Trial_Overlay_Callback;
                end

            end
            hold off;
            if param.dispslices == 1;
                Slice_Display_Callback;
            end
            drawnow expose;
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Refresh statistics
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_statistics

        if size(dim,2)>3 && dim(4)>1 && size(sandbox.hemodynamics,1)==dim(4) ;

            bhemodynamics = (sandbox.hemodynamics+0.2)*(dim(2)/4);
            hemodynamics2 = ([ones(size(sandbox.hemodynamics,1),1) sandbox.hemodynamics])/10;

            hrf_disp = squeeze(idat(Tempvar.slice_l(1),Tempvar.slice_l(2),Tempvar.slice_l(3),:));
            hrf_disp2 = max(hrf_disp)/dim(2);

            if hrf_disp2 >0;
                hrf_disp = hrf_disp/hrf_disp2;
            end
            [best,bint,r,rint,stats] = regress(single(hrf_disp),single(hemodynamics2(:,1:size(hemodynamics2,2))),0.05);
            best=best(2:end);
            bint=bint(2:end,:);
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% refresh overall figure properties
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_figure

        try
            set(0,'CurrentFigure',sandbox_GUI);
        catch
            disp('Main FSB window closed, aborting...')
            return
        end
        dim = size(idat);

        if size(dim,2)>3 && Tempvar.slice_n(4)>dim(4);
            Tempvar.slice_n(4) = 1;
        end
        %set (gcf,'CurrentAxes',txt_win); axis off; axis ij;
        cla (txt_win); cla (beta_win);

        set(h_xaxis,    'Min',1,'Max',dim(1),'Value',Tempvar.slice_n(1),...
            'SliderStep',[(1/dim(1)) (10/dim(1))]);
        set(h_yaxis,'Min',1,'Max',dim(2),'Value',Tempvar.slice_n(2),...
            'SliderStep',[(1/dim(2)) (10/dim(2))]);
        set(h_zaxis,'Min',1,'Max',dim(3),'Value',Tempvar.slice_n(3),...
            'SliderStep',[(1/dim(3)) (10/dim(3))]);

        set(h_xhdr,'String',num2str(sandbox.hdr.dime.pixdim(2)))
        set(h_yhdr,'String',num2str(sandbox.hdr.dime.pixdim(3)))
        set(h_zhdr,'String',num2str(sandbox.hdr.dime.pixdim(4)))

        if size(dim,2)>3 && dim(4)>1;
            set(h_taxis,'Units','normalized','Min',1,'Max',dim(4),'Value',Tempvar.slice_n(4),'SliderStep',[(1/dim(4)) (10/dim(4))]);
            set (hd,'Visible','on');
            if Tempvar.do >1;
                slicet_data = squeeze(idat(:,Tempvar.slice_l(2),Tempvar.slice_l(3),:));
            end
        else
            set (h_taxis,'Visible','off');
            set (hd,'Visible','off');
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% refresh all
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_all


        refresh_data;
        refresh_figure;
        refresh_statistics;
        refresh_image;
        refresh_text;
        refresh_feedb;
        %drawnow expose;

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Refresh statistical map
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_map
        if size(dim,2)> 3 && dim(4)>1 && size(sandbox.hemodynamics,1)==dim(4);
            sandbox.stats.r_cor = zeros(dim(1),dim(2),dim(3),size(sandbox.hemodynamics,2));
            if maskit ==1;
                imean2 = repmat(imean,[1 1 1 size(sandbox.hemodynamics,2)]);
                sandbox.stats.r_cor(imean2<imean_avg/3)=NaN;
            end
            % Mask sandbox.stats.r_cor and all other map matrices with a threshold for brain/nobrain
            r_cor_red = sandbox.stats.r_cor;
            r_cor_neg = sandbox.stats.r_cor;
            sandbox.stats.brob = sandbox.stats.r_cor;
            sandbox.stats.b_est = sandbox.stats.r_cor;
            sandbox.stats.p_rob= sandbox.stats.r_cor;
            sandbox.stats.rob_est= sandbox.stats.r_cor;
            sandbox.stats.rob_prob= sandbox.stats.r_cor;
            sandbox.stats.glm_prob= sandbox.stats.r_cor;
            sandbox.stats.t_test= sandbox.stats.r_cor;
            clustselect = 0;
            % Create a matrix with information on slices already calculated
            Tempvar.iscalc = zeros(dim(1),dim(2),dim(3),15);
            Tempvar.mapcalc = 0;
        end

    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  Refresh volume information
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function refresh_volinfo
        dim = size(idat);
        if size(dim,2)>3
            h = waitbar(0,'Recalculating volume information...');
            hrf_vol_disp = squeeze(mean(mean(mean(idat(1:end,1:end,1:end,:)))));
            waitbar(1/10)
            min_hrf_vol = min(hrf_vol_disp);
            max_hrf_vol = max(hrf_vol_disp);
            avg_hrf_vol = mean(hrf_vol_disp);
            psc_hrf_vol = (max_hrf_vol-min_hrf_vol)/avg_hrf_vol*100;
            Tempvar.ctr_disp2 =0;

            if Tempvar.slice_n(1)>dim(1) || Tempvar.slice_n(2)>dim(2) || Tempvar.slice_n(3)>dim(3) || Tempvar.slice_n(4) > dim(4);
                Tempvar.slice_n = round(dim/2);
                Tempvar.slice_n(4) = 1;
            end
            hdr.dime.dim (2:4) = [dim(1) dim(2) dim(3)];
            waitbar(2/10)
            imean = fsb_calculate_mean(idat);
            waitbar(3/10)
            imean_avg = mean(mean(mean(imean(round(round(dim(1)/2)-round(dim(1)/2)/6):round(round(dim(1)/2)+round(dim(1)/2)/6),...
                round(round(dim(2)/2)-round(dim(2)/2)/6):round(round(dim(2)/2)+round(dim(2)/2)/6),...
                round(round(dim(3)/2)-round(dim(3)/2)/6):round(round(dim(3)/2)+round(dim(3)/2)/6))))); % calculate image average value to use as threshold
            waitbar(5/10)

            psc_vol_string = sprintf('%1.2f',psc_hrf_vol);

            % Intensity based noise area calculation
            % noise_area = single(imean(imean<mean(mean(mean(imean)))/10))

            % Image size based noise area calculation
            noise_lim = round(size(idat)/10);
            noise_lim = noise_lim(1:2);
            noise_area = single(imean(1:noise_lim(1),1:noise_lim(2),:));
            waitbar(6/10)
            %noise_std = std(noise_area(:));
            noise_mean = mean(noise_area(:));
            %Tempvar.snr = imean_avg/noise_std;
            Tempvar.snr = imean_avg/noise_mean;
            waitbar(9/10)
            Tempvar.disp_disc = 0;
            %art_area = idat(dim(1)-0.5*(noise_lim(1):end,:,:,:);
            art_area1 = idat(end-2:end,:,:,:);
            art_mean1 = squeeze(mean(mean(mean(art_area1,3),2),1));
            art_area2 = idat(1:2,:,:,:);
            art_mean2 = squeeze(mean(mean(mean(art_area2,3),2),1));
            art_area3 = idat(round(dim(1)/2),:,:,:);
            art_mean3 = squeeze(mean(mean(mean(art_area3,3),2),1));
            sandbox.art_mean_x = art_mean1+art_mean2;
            sandbox.art_mean_brain = art_mean3;
            %art_area = idat(:,dim(1)-noise_lim(1):end,:,:);
            art_area = idat(:,end-5:end,:,:);
            art_mean = squeeze(mean(mean(mean(art_area,3),2),1));
            %sandbox.discard = zeros(dim(4),1);
            sandbox.art_mean_y = art_mean;
            try
                close(h);
            end
            %set(gcf,'CurrentAxes',he); cla; axis off;
        else
            % This is shitty and needs to be redone
            imean = idat;

        end
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Load parameters to be used
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function reload_params
        dim = size(idat);
        analtrial = 0;
        Tempvar.slice_n = round(dim/2);
        Tempvar.slice_n(4) = 1;
        Tempvar.extr_param = 1;
        %corint_param = 0.1;
        Tempvar.detart_param = 1;
        Tempvar.t_param = 3;
        Tempvar.ssmooth_vox = 2;
        Tempvar.ssmooth_param = Tempvar.ssmooth_vox/2.355;
        Tempvar.mapdisp =0;
        Map.corparam = 0.2;
        Tempvar.do = 1;
        param.displayslices =0;
        %Tempvar.disp_disc = 1;
        idat_spare = idat;
        try
            onset_spare = sandbox.onset;
            intrial_spare = sandbox.intrial;
            hemodynamics_spare = sandbox.hemodynamics;
        end

        best=[]; bint=[];
        hrf_disp =[]; hrf_disp2 = 0;
        bhemodynamics =[];hemodynamics2 =[];
        dctr_mass_xyz =[]; dctr_mass_r =[];

        refresh_buttons;

        if size(dim,2)>3
            if size(dim,2)>3 && dim(4)>1;
                imean = fsb_calculate_mean(idat); % calculate average brain

                imean_avg = mean(mean(mean(imean(round(Tempvar.slice_n(1)-Tempvar.slice_n(1)/6):round(Tempvar.slice_n(1)+Tempvar.slice_n(1)/6),...
                    round(Tempvar.slice_n(2)-Tempvar.slice_n(2)/6):round(Tempvar.slice_n(2)+Tempvar.slice_n(2)/6),...
                    round(Tempvar.slice_n(3)-Tempvar.slice_n(3)/6):round(Tempvar.slice_n(3)+Tempvar.slice_n(3)/6))))); % calculate image average value to use as threshold

                if ~isfield(sandbox,'hemodynamics')

                    if ~exist('hemodynamics','var')
                        sandbox.hemodynamics = [];
                        disp_onset = 0;
                    else
                        sandbox.hemodynamics = hemodynamics;
                        disp_onset = 1;
                    end

                else
                    hemodynamics_spare = sandbox.hemodynamics;
                end

                if exist('ctr_mass','var')==0;
                    ctr_mass = [];
                    Tempvar.ctr_disp2 =0;
                end

                if exist('sandbox','var')==0;
                    sandbox = struct;
                end

                if ~isfield(sandbox,'intrial')
                    if exist('intrial','var')==0;
                        sandbox.intrial = [];
                    else
                        sandbox.intrial = intrial;
                    end
                end
                if ~isempty(sandbox.intrial)
                    if size(sandbox.intrial,2)<3;
                        intrial2 = sandbox.intrial;
                        intrial2(intrial2==0.01)=0;
                        sandbox.intrial(:,3) = bwlabeln(intrial2(:,1));
                    end
                end

                % calculate these variables once to avoid calculating them in the
                % subroutines.
            end

        else
            sandbox.onset =[];
            sandbox.hemodynamics=[];
            hemodynamics_deriv =[];
            intrial =[];
            ctr_mass =[];
            sandbox.discard =[];
            imean = [];

        end

        figure(113); close;
    end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% find and display discard data
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function [sandbox Tempvar] = display_discard(sandbox,Tempvar);
        sandbox.discard = unique(sandbox.discard);
        %sandbox.discard = sandbox.discard;
        trial = sandbox.intrial(sandbox.discard,3);

        trial(trial==0)=[];
        trialtmp =unique(trial);
        trialnum=sandbox.intrial(:,3);
        temp = {};
        temp3=[];

        for xx = 1: size(trialtmp)

            temp{xx} = find(trialnum(:,1)==trialtmp(xx));
            temp2 = cell2mat(temp(xx));
            temp3 = [temp3;temp2];
        end

        temp = temp3;

        sandbox.intrial(sandbox.intrial(:,1)>1) =1;

        if sandbox.intrial(temp,1)==3
            sandbox.intrial(temp,1)=1;
        else
            sandbox.intrial(temp,1)=3;
        end

        disp('----------------------');
        disp(['selected volumes: ' num2str(length(sandbox.discard)) ' of ' num2str(length(sandbox.intrial))...
            ' (' num2str(int8(length(sandbox.discard)/length(sandbox.intrial)*100)) '%)']);
        disp('--');
        disp(num2str(sandbox.discard'));
        disp('----------------------')
        disp(['selected trials: ' num2str(length(trialtmp)) ' of ' num2str(max(sandbox.intrial(:,3)))...
            ' (' num2str(int8(length(trialtmp)/max(sandbox.intrial(:,3))*100)) '%)']);
        disp('--');
        disp(num2str(trialtmp'));
        disp('----------------------');
        disp(['Trial / volume ratio : ' num2str(length(trialtmp)/length(sandbox.discard))]);

        if ~isempty(sandbox.discard);
            set(gcf,'CurrentAxes',hd);
            hold on; plot (sandbox.discard,10,'-w+',...
                'LineWidth',1,...
                'MarkerEdgeColor','w',...
                'MarkerFaceColor',[1 1 1],...
                'MarkerSize',10); hold off;
            Tempvar.disp_disc = 1;

            Tempvar.feedb =['Artifact detection threshold ' num2str(Tempvar.detart_param) ];

        else
            disp('Error detecting artefacts');
            Tempvar.feedb = ('Error detecting artefacts');
            refresh_feedb;
            Tempvar.disp_disc =0;
            sandbox.discard = [];
        end


    end

end



